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ABSTRACT 


Electromagnetic  scattering  from  trees  and  vegetation  is  of  prime  importance  in 
radar  and  remote  sensing.  The  actual  problem  of  scattering  from  trees  is  rather 
complicated  and  involves  three  dimensional  scattering  from  lossy,  electrically  large, 
and  randomly  oriented  objects. 

In  this  thesis,  the  radar  cross  section  of  a  planar  fractal  tree  is  considered. 
Although  a  planar  tree  is  far  from  being  real,  scattering  from  it  sheds  light  on  the 
scattering  phenomenon  from  an  actual  tree.  The  planar  tree  is  generated  using 
fractal  geometry  and  its  branches  are  considered  perfectly  conducting.  The  tree  is 
illuminated  by  a  plane  wave  and  the  problem  is  solved  using  the  moment  method. 
Data  is  presented  for  the  radar  cross  section  for  different  branching  angles  of  the 
tree  and  at  different  frequencies. 
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I.  INTRODUCTION 


A.  NEED  FOR  THE  STUDY 

The  trees  existing  in  the  natural  world  are  fractal  anisotropic.  They  are  made 
up  of  long,  intersecting  and  lossy  objects.  The  geometry  of  these  objects  is  not  easy 
to  set  up  as  they  are  randomly  oriented  in  a  three  dimensional  space.  The  use  of 
fractals  facilitates  the  modeling  of  semi-randomly  distributed  structures. 
Mandelbrodt  [Ref.  1]  has  shown  that  a  number  of  naturally  occurring  phenomenon 
such  as  coastlines,  clouds,  trees,  etc.  are  fractal  in  nature.  For  instance,  when  a 
branch  is  divided  into  two  (or  more),  the  ratio  of  the  length  of  the  subbranches  to 
the  main  branch  length  remains  constant.  Furthermore,  the  branching  angle  also 
remains  the  same. 

In  this  thesis,  scattering  from  planar  fractal  trees  is  considered.  The  radar 
cross  section  of  a  real  fractal  tree  is  a  complicated  scattering  problem.  The  analysis 
of  this  problem  in  a  two  dimensional  space  gives  an  approximate  idea  of  the  radar 
cross  section  of  a  real  tree. 

The  radar  cross  section  of  au  object  is  a  quantitative  measure  of  the  ratio  of 
the  power  density  that  is  received  and  scattered  by  the  object  to  the  power  density 
of  the  electromagnetic  wave  that  illuminates  that  object.  The  radar  cross  section  is 
independent  of  the  range  of  the  object  for  the  far— field  situation.  The  theoretical 
definition  of  the  radar  cross  section  "ou  is  given  by  the  formula: 

<>  |  | 

a  =  4*R"  Urn  — \ 

R-  a.  1  £'  i 
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where  fi1  is  the  incident  electric  field  vector,  fis  is  the  scattered  electric  field  vector, 
and  R  is  the  distance  between  the  scattered  object  and  the  point  of  observation.  The 
radar  cross  section  has  dimensions  of  area.  Usually,  it  is  expressed  in  square 
wavelengths. 

B.  STATEMENT  OF  THE  PROBLEM 

This  thesis  investigates  the  radar  cross  section  of  planar  fractal  trees.  These 
trees  are  composed  of  planar  thin  strip  dipoles  of  arbitrarily  dimensions  and 
orientations  in  a  plane.  These  structures  are  excited  by  plane  waves  of  various 
frequencies. 

The  first  step  in  solving  the  problem  is  to  calculate  the  induced  current 
distribution  on  each  strip.  The  calculation  of  the  current  distribution  is  based  on  the 
theory  of  the  moment  method  and  requires  a  knowledge  of  the  impedance  between 
any  two  of  these  strips  as  well  as  the  voltage  on  each  planar  strip  due  to  the 
incident  electric  field.  The  basic  concepts  and  the  calculation  of  the  current 
distribution  are  described  in  Chapter  2. 

In  Chapter  3,  the  development  of  a  FORTRAN  program,  is  presented.  The 
evaluation  of  the  radar  cross  section  requires  the  knowledge  of  the  scattered  electric 
field  due  to  the  induced  currents  on  the  planar  strips.  The  program  computes  the 
scattered  electric  field  and  then  the?  radar  cross  section  of  that  structure.  The  details 
of  these  calculations  are  also  presented  in  this  Chapter. 

The  computer  models  that  are  used  by  the  develojred  program  are  presented  in 
Chapter  4.  Their  generation  is  based  on  the  fractal  geometry.  An  existing  and 
modified  program  is  used  to  generate  the  geometry  of  the  planar  fractal  trees  in 
order  to  be  used  as  input  in  the  developed  program. 


The  numerical  results  of  the  radar  cross  section  of  a  single  planar  dipole  and  a 
a  number  of  planar  fractal  trees  are  presented  in  Chapter  5.  The  scattering  from  a 
single  dipole  is  compared  with  standard  results  for  a  similar  case.  The  limitations  of 
the  developed  program  are  also  presented  in  this  Chapter. 

In  Chapter  5,  the  conclusions  of  the  radar  cross  section  results  and 
recommendations  are  presented.  The  two  programs  that  are  used  for  this 
investigation  are  listed  in  the  Appendices. 
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II.  METHOD  OP  MOMENTS  THEOItY 


In  tills  section  of  the  study,  the  basic  concepts  of  the  moment  method  theory 
are  presented.  This  theory  is  used  in  the  development  of  the  ItCS  program  to  find 
the  current  distribution  on  a  planar  strip  due  to  an  incident  plane  wave. 

For  a  given  structure  consisting  of  planar  dipoles  the  impedance  between  any 
two  of  them  is  calculated  from  the  knowledge  of  the  geometry  and  the  wavelength  of 
the  incident  plane  wave.  The  voltage  on  each  dipole  is  calculated  from  the 
knowledge  of  the  characteristics  of  the  incident  electric  field.  The  induced  current 
distribution  on  each  dipole  of  the  structure  is  determined  from  the  calculated 
impedance  and  voltage  using  the  method  of  moments  theory. 

A.  GENERAL  THEORY 

The  method  of  moments  is  a  numerical  procedure  for  solving  integral 
equations  of  the  form: 

b 

f(x')K(x,x')dx  at  g{x),  a  <  x  <  b  (eqn  2.1) 

••a 

where  f(x')  is  an  unknown  function,  K(x,x')  is  a  known  Kernel  or  Green's  function, 
and  g{x)  is  a  given  function.  This  procedure  reduces  the  integral  equation  (eqn  2.1) 
to  a  system  of  simultaneous  linear  algebraic  equations  in  terms  of  some  unknown 
coefficients.  This  method  requires  that  the  function  f(x')  be  approximated  by  a 
series  of  N  expansion  functions  or  "  basis  functions  ",  such  that 


n=l,2, . ,N 


(eqn  2.2) 


f(x')  ~  £  anfn(x'), 

n«l 

where  the  domain  of  fn(x')  is  the  same  as  that  of  f(x')  and  an's  are  the  complex 
unknown  expansion  coefficients. 

There  are  two  types  of  basis  functions.  The  subdomain  functions,  which  are 
nonzero  over  a  part  of  the  domain  of  the  unknown  function  f(x'),  and  entire  domain 
functions  being  nonzero  over  the  entire  domain  of  f(x').  In  antennas,  some 
commonly  employed  subdomain  basis  functions  are  the  piecewise  sinusoid  functions, 
the  unit  height-pulse  functions  and  the  piecewise  triangular  functions. 

The  subdomain  procedure  requires  subdivision  of  the  structure  into  N 
nonoverlapping  segments.  Figure  2.1  shows  a  segmented  line  where  the  segments  are 
assumed  to  be  collinear  and  of  equal  length,  although  this  condition  is  not  necessary. 


X(  x^  xj  x*  x$  b 


Figure  2,i  Segmented  Line.  [From  Ref.  2| 
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Figure  2.2  shows  a  subdomain  unit  height-pulse  function  which  produces  a 
staircase  representation  of  the  unknown  function  f(x').  Figure  2.3  shows  a 
subdomain  sinusoid  basis  function  and  the  representation  of  the  function  f(x'). 
Figure  2.4  shows  a  subdomain  triangular  basis  function  producing  a  smoother 
representation  of  the  function  f(x')  than  the  case  of  the  unit  height-pulse  basis 
function. 

The  use  of  entire-domain  basis  functions  does  not  require  any  segmentation  of 
the  structure.  One  of  the  most  most  commonly  used  basis  functions  of  this  kind  is 
the  sinusoidal  basis  functions. 
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X 


Figure  2.2  Unit  Height-Pulse  Basis  Function  and  a  Staircase 
Representation  of  f(x').  [From  Ref.  2] 
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Figure  2.3  Sinusoid  Basis  Function  and  a  Representation 
of  fix').  (From  Ref.  3] 


Figure  2.4  Triangular  Basis  Function  and  a  Representation 
of  f(x').  [From  Ref.  21 
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The  substitution  of  the  function  f(x')  by  a  sequence  of  N  basis  functions  leads 
to  one  equation  of  N  unknowns  of  the  expansion  coefficients  an,  which  can  be  found 
by  using  N  linearly  independent  equations.  These  equations  are  set  up  by  the  use  of 
"testing  or  weighting"  functions  <Jm(x),  m=l,2,....N.  At  this  point  the  definition  of 
the  inner  product  is  required.  The  inner  product  (f,  g)  between  any  two  functions  or 
vectors  f,  g  is  a  scalar  operation  defined  as 


<i,g)=  f'gds 


(eqn  2.3) 


where  S  is  the  surface  of  the  structure  that  is  analyzed  [Ref  4].  The  inner  product  of 
the  selected  testing  functions  u4(x),  with  the  two  sides  of  the  original  integral 
equation  leads  to  the  equation: 


r(x/)K(x,x/)dx/  ^  -  ^g(x),  <*(x)  ^ 


(eqn  2.4) 


or 


(  ufc(x),  l  a„f.>(x')dx'  y  =  g(x),  u4(x)  ^ ,  m  -  1,2, . ,N 

J  a  n= 1 

(eqn  2.5) 

The  use  of  the  above  definition  for  the  inner  product,  yields: 


N  r  b  j*  b  r  b 

£a„  ^(xjdx  f„(x')K(x,x')dx'  =  g(x)oim(x)dx,  m  =  1,2 . N 

„=i  Ja  Ja  ja 

(eqn  2.6) 
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(eqn  2,7) 


The  following  substitutions 

T  b 

VB  =  g(x)u;m(x)dx,  m=l,2, . ,N 

Ja 

and 


Z 


mn  — 


'  b 

a4(x)dx 

Ja 


b 

fn(x/)K(x,x/)dx/, 

Ja 


m,n=l,2,....N 


result  in  a  matrix  equation  of  the  form: 


(eqn  2.8) 


[Zmn]  *  fan]  —  (V m]  (oqn  2.9) 

The  unknowns  [an]  can  be  obtained  through  a  matrix  inversion: 

fan]  =  l^mn]  '  [Vm]  (eqn  2.10) 

where  [Zmnj  is  the  inverse  matrix.  The  column  vector  [Vm]  depends  upon  the  given 
function  g(x)  and  the  selected  testing  functions  ^(x).  The  matrix  [Zmn]  depends 
upon  the  known  kernel  K(x,x')  and  both  the  selected  basis  and  testing  functions. 
Once  the  expansion  coefficients  are  known,  the  function  f(x')  is  also  known. 

The  choice  of  basis  and  testing  or  weighting  functions  is  based  upon  experience 
and  the  rule  is  that  their  number  has  to  be  the  same.  The  procedure  of  using  the 
same  basis  and  testing  functions  is  called  the  "  Garlekin's  method  ". 


B.  APPLICATIONS  TO  EM  THEORY 

A  number  of  problems  in  electromagnetic  radiation  and  scattering  can  be 
solved  by  the  method  of  moments.  In  the  present  case,  the  simple  case  of  a  perfectly 
conducting  object  "  situated  in  free  space  "  is  considered. 

The  basic  problem  is  to  investigate  the  case  when  an  object  is  illuminated  by 
fields  of  known  impressed  electric  and  magnetic  currents  (J1,  M1).  In  the  absence  of 
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the  object,  the  impressed  currents  radiate  the  assumed  known  incident  electric  and 
*  • 

magnetic  fields  (E1,  H1).  In  the  presence  of  this  object,  the  impressed  currents 
radiate  the  unknown  total  fields  (El,  H*). 

The  integral  equation  is  obtained  by  the  surface  equivalence  principle, 
replacing  the  object  by  free  space  together  with  the  electric  surface  current  density 

J  =  nxHt  (eqn2.11) 

where  J  exists  on  the  entire  surface  S  of  the  object  and  n  is  the  unit  vector  normal 
to  that  surface.  In  free  space,  J  radiates  the  scattered  fields  Es,  Hs 

E8  =  E*  —  E*  (eqn  2.12) 

H8=Ht->Hi  (eqn  2.13) 

The  boundary  conditions  for  this  case  enforce  the  total  tangential  electric  field 
on  the  surface  S  to  zero: 

nx(Es  +  Ei)  =  0  (eqn  2.14) 

This  is  an  integral  equation  for  J  since  the  scattered  electric  field  Es  can  be  written 
as  an  integral  over  S  of  the  dot  product  of  J  and  the  dyadic  free  space  Green's 
function.  (Ref.  5] 

As  the  geometry  of  the  object  is  known,  it  is  more  convenient  to  use  the  total 
current  I  instead  of  the  total  current  density  J.  So,  the  problem  is  to  find  the 
unknown  current  I  induced  by  the  incident  electric  field.  The  method  of  moments 
solve  these  kind  of  problems  by  the  following  procedure. 

The  first  step  is  to  expand  the  unknown  current  I  in  terms  of  some  basis  set: 

N 

la  SIaFu  (eqn  2.15) 

ns| 

where  the  lu  are  the  sequence  of  N  unknown  complex  coefficients,  and  the  Fu  is  a 
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sequence  of  N  known  modes  or  basis  functions.  The  best  choice  of  Fn  for  a  given 
problem  could  be  quite  involved  and  is  discussed  in  [Ref.  4,  pp  308-310]. 

The  second  step  is  to  select  the  testing  or  weighting  functions  m=l,2,....N. 
These  can  be  identical  with  the  basis  functions  or  different.  The  inner  product  of  the 
sequence  of  N  weighting  functions  ufo,  with  both  sides  of  the  integral  equation  gives 
a  N  x  N  system  of  simultaneous  linear  algebraic  equations  of  the  symbolic  form: 

[Z].[I]  =  [V]  (eqn  2.16) 

where  I  is  the  current  column  vector  whose  N  components  give  the  values  of  ln 
and  [Z]  is  the  N  x  N  impedance  matrix  given  by  the  equation 


(eqn  2.17) 


where  S  is  the  surface  of  the  structure  being  analyzed.  The  impedance  matrix  [Z]  is 
always  symmetric,  and,  for  the  special  case  of  a  thin  dipole  instead  of  an  arbitrary 
object,  is  also  a  Toeplitz  matrix.  In  a  Toeplitz  matrix,  Zmn  depends  only  on  |m-n  | . 

Generally,  [Z]  is  dependent  only  on  the  geometry  and  material  composition  of 
the  scatterer,  but  not  on  the  incident  fields  [Ref  5).  The  right-hand  side  of  the  last 
equation  is  the  voltage  vector  whose  N  components  give  the  corresponding  mode 
voltage.  The  voltage  vector  depends  only  on  the  excitation,  i.e.,  the  incident  electric 
field.  The  dimensions  of  the  elements  of  [Z]  and  [V]  are  volt-amps  (VA),  while  the 
elements  of  [I]  are  dimensionless. 

The  solution  of  the  last  matrix  equation  is  the  column  vector  [1],  whose 
,'lements  represent  the  complex  coefficients  I0.  As  the  total  current  I  was  expanded 
by  a  sequence  of  N  known  modes  or  basis  functions  and  the  complex  coefficients  1„ 
are  known,  the  total  current  and  so  the  total  current  density  J  is  also  known. 
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Although  the  choice  of  weighting  functions  is  free,  it  has  to  be  considered  that 
the  matrix  equation  being  solved  requires  the  evaluation  of  N  terms  and  each  term 
requires  two  or  more  integrations.  When  these  integrations  are  to  be  done 
numerically,  the  computations  become  complicated.  There  is  a  way  to  reduce  this 
complexity  by  choosing  as  weighting  functions  the  ihrac  delta  functions.  This  is  the 
method  of  point-matching  in  which  delta  functions  are  enforced  only  at  discrete 
points  on  the  surface  S.  The  results  of  this  method  can  be  quite  accurate  especially 
when  the  discrete  points  are  selected  to  be  equally  spaced.  The  solution  satisfies  the 
electromagnetic  boundary  conditions  (e.g.,  vanishing  tangential  electric  fields  on  the 
surface  of  an  electric  conductor)  only  at  discrete  points  [Ref  4].  Between  these 
points  the  boundary  conditions  may  not  be  satisfied.  In  this  case  it  is  required  to 
define  the  deviation  as  a  residual  (e.g.,  residual=AE|tan  =  E  |tan+  E  [tan  t  0  on 
the  surface  S  of  an  electric  conductor)  and  use  the  method  of  weighted  residuals  so 
that  the  boundary  conditions  will  be  satisfied  in  an  average  sense  over  the  entire 
surface  S. 
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III.  ANALYSIS  AND  DEVELOPMENT  OF  RCS  PROGRAM 


In  this  chapter,  the  basic  steps  that  are  involved  in  the  calculation  of  the 
Radar  Cross  Section  (RCS)  of  a  planar  structure  composed  of  conducting  strips  are 
presented.  As  mentioned  earlier,  the  fractal  tree  is  modeled  as  consisting  of  planar 
thin  strips.  The  radar  cross  section  of  the  fractal  tree  is  calculated  using  the  method 
of  moments.  A  Fortran  program  is  developed  to  calculate  the  RCS  of  the  tree  for  a 
specified  geometry  and  at  a  given  frequency. 

The  moment  method  discretizes  an  integral  equation  to  a  matrix  equation  of 
the  form: 

|Z]  •  [I]  =  [V]  (eqn  3.1) 

where  [Z]  is  the  impedance  matrix  whose  elements  represent  the  mutual  impedance 
between  any  two  dipoles  of  a  given  structure  depending  upon  the  wavelength  and 
geometry  of  the  structure,  (Vj  is  the  voltage  matrix  whose  elements  correspond  to 
the  voltage  on  each  dipole  due  to  excitation  ( incident  electric  field  ),  and  (Ij  is  the 
unknown  matrix  whose  elements  represent  the  induced  current  on  each  dipole. 

A.  DEVELOPMENT  OF  RCS  PROGRAM 

The  numerical  results  for  RCS  are  obtained  from  a  Fortran  program  that 
calculates  backscattered  RCS  of  a  planar  structure  consisting  of  arbitrarily  oriented 
dipoles.  For  convenience,  the  structure  will  be  assumed  to  be  in  the  y-z  plane  of  an 
xyz  cartesian  coordinate  system. 

Figure  3.1  shows  the  structure  to  be  investigated.  A  large  number  of  planar 
dipoles  of  variable  lengths,  widths,  and  orientations  in  the  y-z  plane  is 


15 


illuminated  by  a  plane  wave  with  electric  field  linearly  polarized  and  characterized 
by  the  angles  <po  and  6q.  It  is  required  to  calculate  the  backscattered  RCS  from  this 
structure. 


The  incident  electric  field  £*  is  given  by  the  formula 


£i  _  0.e— j(kx*x  +  ky*y  +  k2*z) 

(eqn  3.2) 

where  Ic  is  the  propagation  vector  and 

kx  +  ky  +  k2  =  ko  =  J /io£o 

(eqn  3.3) 

k0=  2jr 

A 

(eqn  3.4) 

kx  =  kosin^ocosyjo 

(eqn  3.5) 

ky  =  kosinflosinv*) 

(eqn  3.6) 

kz  =  kocostfo 

(eqn  3.7) 

e  =  Exx  +  Eyy  +  Ezz 

(eqn  3.8) 

The  electric  field  fi1  lies  on  the  plane  perpendicular  to 

the  direction  of 

propagation  of  the  plane  wave.  Hence,  e*k  =  0  implies  that 

Exkx  +  Eyky  +  E2kz  —  0 

(eqn  3.9) 

The  substitution  of  Ey  and  E2  by  the  variables  a  and  b  (independent  variables)  gives 


the  coordinate  Ex  as  a  dependent  variable 


(akv  -4-  bka) 

kx 


(eqn  3.10) 


Equation  3.1  the  becomes 


ay  +  bz  -  t  aky  *  )  x  c~^  +  +  ) 

kx 


(eqn  3.11) 
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The  RCS  program  that  has  been  developed  makes  the  following  assumptions: 

1.  The  dipoles  are  very  thin  planar  strips.  The  width  of  the  dipoles  is 
assumed  to  be  electrically  small  so  that  only  the  axial  current  is 
significant.  Further,  there  is  no  current  variation  along  the  width  of 
the  dipoles. 

2.  The  dipoles  are  not  intersecting. 

3.  The  dipoles  are  perfect  conductors. 

The  program  uses  the  theory  of  moment  method  and  the  way  that  it  solves  the 
problem  can  be  understood  if  a  single  planar  thin  dipole  in  the  y-z  plane  is 
considered.  Each  dipole  is  subdivided  into  equal  segments.  Overlapping  Piecewise 
Sinusoidal  (PWS)  modes  are  assumed  to  exist  on  the  dipole.  The  length  of  each 
PWS  mode  is  equal  to  the  length  of  two  segments.  Figure  3.2  shows  the  mth  PWS 
mode.  The  PWS  mode  has  a  length  2h0,  and  a  width  2wm.  The  coordinates  of  its 
center  are  (yB,  zm),  and  it  is  oriented  at  an  angle  measured  from  the  i  axis. 

The  incident  electric  field  induces  current  along  the  aids  (  of  that  mode.  In 
Figure  3.2,  the  induced  current  J  varies  sinusoidally  along  the  axis  of  the  mu,  mode. 
A  new  coordinate  system  is  introduced,  where  ( is  the  axis  of  the  dipole  and  7/  is 
the  axis  perpendicular  to  (.  The  coordinate  system  is  obtained  by  rotating  the 
y-z  system  about  the  x-axis  through  an  angle  90°-t%.  Figure  3.3  shows  the  details 
of  this  rotation. 

The  relation  between  the  coordinates  of  the  center  of  the  mode  along  the  axis 
of  the  two  orthogonal  systems  is  found  to  be: 

>*  -  .Vo  +  Csin(&,)  -  qcos(ifa)  (oqn  3.12) 

z  =  za  +  (cos(&.)  +  i$in(tfe)  (eqn  3.13) 


IS 


>> 

♦ 


O 


Figure  3.3  Traosfornuitio:.  of  PW$  Mode's  Center  Coordinates. 
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1.  Voltage  Equations 

In  this  section,  expressions  for  the  elements  of  the  voltage  matrix  [V] 
due  to  the  incident  electric  field,  are  presented.  Each  element  of  this  matrix 
corresponds  to  a  PWS  mode  of  the  structure  that  is  investigated.  The  size  of  the 
voltage  matrix  is  equal  to  the  total  number  of  modes  of  the  structure. 

The  induced  current  density  on  the  mth  mode  of  the  structure,  due  to 
incident  electric  field,  is 


J  =  J-  (eqn  3.16) 

2wm  sin(k0hm) 

where  2wm  is  the  width  and  2hm  is  the  length  of  the  mth  mode.  For  each  mode,  the 
induced  current  variation  will  be  sinusoidal  along  its  axis,  being  maximum  at  the 
center  and  zero  at  the  end  points.  The  corresponding  voltage  Vm  is 


'  u>m  r  hm 
— ~ hm 


£’  •  J  d(d T) 

!  x=o 


(eqn  3.17) 


From  equations  3.11  and  3.16,  it  is  seen  that 


E1 


•3 


[a  sin)*,)  +  b  cos(*)|  e"^*  +  +  Mz"  +Ccos(^„)l 


(eqn  3.18) 

where  ipm  is  the  angle  between  the  z  axis  (reference  for  angle  measurements)  and  the 
axis  of  the  mode  and  yra,  zK  are  the  coordinates  of  the  center  of  the  mode  along  the 


21 


y  and  z  axis  respectively.  The  substitution  into  equation  3.17,  gives  the  following 
form  of  Vm: 


Vm  = 


p-j(kyym  +  kzzm) 


hm 


r 

[a  sin(fe)  +  b  cos(*,)]  J  e-*;<  sin(kt,(h„-|  (|  ))d< 


sin(kohm)  ’  ”  J-hm 

(eqn  3.19) 

where  k<-  =  kysin(^ra)  +  kzcos(^>m).  A  closed  form  evaluation  of  the  integral  in  this 
equation  leads  to  final  form  of  Vm 


\r  _  2k0Vora 

Vm  —  — ^ — ^—7r 

kr  -  ko 


cos(k0hm)  -  cos(kf  hm)  j 


Vm  —  Vom  hro  sin(kollm) 


where 


k.Q  ^  i  ko 

(eqn  3.20) 
k^  ~  ^  ko 
(eqn  3.21) 


Vom  — 


e—j(kyyra  +  kzzm) 
sin(kohm) 


a  sin($m)  +  b  cos(0m) 


(eqn  3.22) 


2.  Radiation  Equations 

In  this  section,  expressions  for  the  far-zone  scattered  fields  due  to  the 
mo,  PWS  mode  are  developed.  For  far— field  observations  the  electric  field  is  given 
in  spherical  coordinates  by  the  following  equations  (Ref  4): 


(eqn  3.23) 
(eqn  3.24) 


Ef  a  0 


E 


■kne~jk»r 


(L*>+”N«) 


4m: 


where 


(eqn  3.25) 


L«=.. 

[Mxcosflcos^>  +  Mycosfein^-Mzsintf]  e+^°r  cos^“ 
S 

ds' 

(eqn  3.26) 

% 

[-Mxsin^  +  Mycosd  e+^r  cos^m  ds' 

JS 

(eqn  3.27) 

>v.' 

[Jxcos0cos<^  +  Jycosfeiny?  -  J2sin0]  e^01  cos^m  ds' 

JS 

(eqn  3.28) 

vJ 

[— JxSinv?  +  Jycostjp]  e^o1  cos^m  ds' 

s 

(eqn  3.29) 

Tj  =  120* 

(eqn  3.30) 

The  quantities  Jx,  Jy,  Jz,  are  the  components  of  the  electric  current 
density  3S  that  are  induced  on  the  mtb  mode  over  the  surface  S,  and  the  quantities 
Ms,  My,  M2>  are  the  coordinates  of  the  magnetic  current  density  over  the  surface 
S.  As  the  structure  is  on  the  y-z  plane  and  $s  is  zero,  the  quantities  Jx,  L^,  are 
zero.  Equations  3.24  and  3.23  become: 


E 


-ihcc-jk»r 
4  <rr 


4  *T 


% 

(eqn  3.31) 

(eqn  3.32) 
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As  the  current  density  is  along  the  (-axis,  all  quantities  within  the 
integrals  for  3  and  K4  will  be  expressed  in  terms  of  the  rj-£  system.  The  surface 
element  that  is  used  in  all  integrals  is  ds'  =  dydz  =  d(dr?  =  d(,  as  the  width  of  the 
mode  is  assumed  to  be  very  small  in  terms  of  its  length.  The  transformations  from 
the  y-z  system  to  the  (-77  system  are  made  by  using  the  following  substitutions: 

r/cos(^m)  =  ysinfcin <p  +  zcos0  (eqn  3.33) 

where 

y  =  ym  +  (sin(^ni)  (eqn  3.34) 

z  =  zm  +  (cos(V'm)  (eqn  3.35) 


ana 


—  dyy  "h  JzZ  —  ( 


In 


2wmsin(kohm) 


sin(k0(hm-|(|)) 


where 


(eqn  3.36) 


Jy  =  J^sin(^m) 

Jz  =  J  ^  cos(^m) 

These  substitutions  and  algebraic  manipulations  give  the  quantities 
for  the  mth  mode  in  the  (-77  system: 


(eqn  3.37) 

(eqn  3.38) 

N«  and  N 
Um  <^m 


N 


Om 


-2k  0 
1 - 1 


-  ko 


cos(Emhni)  -  cos(kohm) 


N0m  =  Nfl0  h°  sin(k°h») 


Em  f  iko 


(eqn  3.39) 

Em  —  ±ko 

(eqn  3.40) 


and 
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Ejn  ±ko 


-2k  o 

2 - 7 

Ea  .  k0 


COS(Enjllin)  -  COS(kohm) 


N(pm  =  h»  sin(Effihm) 

where 


(eqn  3.41) 

Em  = 

(eqn  3.42) 


^  _  sin^mcos feint/?  -  cos^msinfl  ^  ^[ko^mSin&ini/H-ZmCOS#)] 

sin(kohm) 


N  _  sin^mcos(p  j  ej[ko(ymsinfein^+zmcos^)] 
sin(kohm) 


(eqn  3.43) 


(eqn  3.44) 


Em  =  ko(sin^n,sin$3in^  +  cos^mcosfl) 

(eqn  3.45) 

In  this  thesis,  only  the  monostatic  radar  cross  section  will  be 
considered.  In  the  radiation  equations  the  angles  6  and  y  represent  the  orientation 
angles  of  the  scattered  electric  field  £s.  Their  relation  with  the  incident  angles  Oq 
and  tpo  for  the  monostatic  case  is: 

0  =  7r  -  0o  (eqn  3.46) 

V?  =  7T  —  c^o  (eqn  3.47) 

If  M  denotes  the  total  number  of  PWS  modes  in  the  structure  under 
investigation,  then 

N/i  _  E  (eqn  3.4S) 

u  ~  m=l  wn 


N  =  £  N  (eqn  3.49) 

*  in=l  • 
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The  final  expression  for  the  equations  3.24  and  3.25  is 


E«=c-  s,Nto 
m=l 

u 

E  =  C  S  N 
*  m=l  ^ 


where 


4xr 


(eqn  3.50) 
(eqn  3.51) 


(eqn  3.52) 


3.  RCS  Equations 

When  an  incident  plane  wave  with  electric  field  fi1  strikes  the  object 
and  fis  is  the  scattered  electric  field,  the  radar  cross  section  o  is  defined  as 


a  =  lim  4  7rr 
r  -*  oo 


,2  |fisl 


(eqn  3.53) 


For  the  scattered  electric  field  £s,  its  spherical  coordinates  E^,  and  E^  have  been 
calculated  by  the  equations  3.50  and  3.51.  The  incident  electric  field  £'  is  known  by 
means  of  equation  3.11. 


2  2 
a|  +  |bj  + 

aky  bkj 

2 

kx 

2  2 

2  r 

2  2-, 

+|E„I 

=  iC|  [[|N,|  +|NV|  j 

(eqn  3.54) 


(eqn  3.55) 
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The  combination  of  equations  3.11,  3.50-3.52,  3.54-3.55,  and  3.53  leads  to  the  final 
formula  for  RCS: 


4?r 


|n/+|n/ 


1 


|a|  +  |b|  4 


'  4  bkz 


(eqn  3.56) 


This  is  the  formula  that  the  program  uses  to  compute  the  RCS  of  a 
planar  structure  for  a  given  set  of  incident  angles  $o,  and  y?o. 


B.  INPUT-OUTPUT  OF  RCS  PROGRAM 

The  RCS  program,  which  is  described  in  Appendix  A,  is  a  FORTRAN 
program  having  two  input  files.  The  first  input  file,  INPUT  1 ,  contains  the  data  that 
characterize  the  incident  plane  wave  and  the  data  that  describe  the  geometry  of  the 
planar  structure  whose  broadside  RCS  is  measured.  The  second  input  file,  INPUT, 
contains  the  set  of  incident  angles  0o,  <A)- 

The  program  reads  from  the  INPUT1  file  the  following  input  data: 

1.  frequency  of  the  incident  plane  wave  in  GHz, 

2.  parameters  a  and  b  that  characterize  the  polarization  of  the  incident 
electric  field. 

3.  number  of  dipoles  that  the  structure  consists  of, 

4.  half  length  of  each  dipole  in  cm, 

5.  half  width  of  each  dipole  in  cm, 

6.  coordinates  in  cm  of  the  center  of  each  dipole  along  the  two  axes  of 
the  orthogonal  system  that  the  structure  lies, 
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7.  orientation  angle  of  each  dipole,  measured  from  the  vertical  axis, 
positive  in  the  clockwise  direction,  and 

8.  number  of  segments  of  each  dipole. 

As  the  program  reads  these  input  data  it  generates  the  geometry  of  the  PWS 
modes,  calculates  the  impedance  elements  between  any  two  modes  of  the  structure, 
and  fills  the  impedance  matrix  [Z].  The  size  of  the  impedance  matrix  depends  on  the 
total  number  of  modes.  When  this  matrix  is  calculated  and  filled,  the  program 
inverts  it  to  [Z]_1  and  stores  it  for  later  use. 

The  next  step  is  to  read  from  the  INPUT  file  the  sets  of  incident  angles  6q,  <po 
and  for  each  one  it  calculates  the  voltage  on  each  mode,  filling  the  voltage  matrix 
[V].  This  matrix  is  a  column  vector  which  depends  upon  the  excitation  only.  Then  it 
multiplies  the  stored  inverted  impedance  matrix  [Z]-1  by  the  voltage  matrix.  This 
yields  the  current  vector  [I]: 

[I]  =  [zrMv]  (eqn  3.57) 

As  the  induced  currents  are  known,  the  program  uses  the  previously  described 

radiation  equations  and  calculates  the  RCS  of  the  structure  corresponding  to  the 

given  set  of  incident  angles  Oo ,  <A).  The  output  RCS  is  normalized  to  the  square  of 
o 

the  wavelength  A  ,  and  is  given  in  dB.  The  program  reads  the  next  set  of  incident 
angles  6q ,  <pa  and  repeats  the  same  procedure  to  compute  the  RCS  of  the  new  set. 

Although  the  input  lengths  and  widths  are  in  cm,  the  program  considers  them 
normalized  to  the  wavelength.  For  accurate  results  at  least  4  segments  per 
wavelength  are  chosen.  The  selection  of  the  width  of  each  dipole  is  arbitrarily  taken 
as  L/W  =  33. 


2S 


IV.  FRACTAL  TREE  GENERATION 


The  problem  of  measuring  the  radar  cross  section  of  a  real  natural  tree  is 
complicated  as  it  requires  an  investigation  in  three  dimensional  space  with  very 
long,  intersecting  elongated  objects  that  are  randomly  oriented. 

For  simplicity  the  radar  cross  section  (RCS)  problem  is  solved  in  a  two 
dimensional  space.  RCS  is  calculated  from  a  planar  fractal  tree  whose  branches  are 
considered  as  thin  planar  dipoles  of  different  lengths  and  widths.  In  an  actual  tree, 
the  branches  are  lossy  and  in  general  anisotropic.  However,  in  the  present  model, 
the  tree  is  comprised  of  perfectly  conducting  branches.  The  geometry  of  this  tree  is 
based  upon  the  fractal  geometry. 

A.  DEFINITIONS 

Fractal  is  a  mathematical  set  or  object  whose  form  is  extremely  irregular  or 
fragmented  at  all  scales  [Ref.  6].  The  requirement  to  describe  the  shape  of  many 
objects  that  appear  in  the  natural  world,  such  as  trees,  mountains,  coastlines,  etc., 
led  to  the  generation  of  fractal  geometry.  As  Euclidean  geometry  cannot  give 
mathematical  expressions  to  describe  fractal  objects,  fractal  geometry  is  used  to 
describe  mathematically  many  natural  patterns. 
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The  basic  characteristics  of  fractal  objects  are  [Ref.  7]: 

1.  A  large  degree  of  heterogeneity. 

2.  A  self-similar  structure  over  many  size  scales.  Self-similarity  refers 
to  the  general  preservation  of  form  or  characteristic  regardless  of  the 
scale  of  observation. 

3.  The  lack  of  a  well-defined  (characteristic)  scale. 

The  geometric  characteristics  of  fractal  objects  are  useful  for  describing  phenomena 
of  nature,  such  as  scattering  of  objects  like  landscapes  or  surface  cracks.  Although 
these  objects  are  irregular  and  randomly  oriented  in  nature,  they  show  structural 
similarities  on  several  different  discrete  size  scales  [Ref.  7]. 

One  measure  of  structural  complexity  is  the  fractal  dimension  Df.  There  are 
several  definitions  of  Df  depending  upon  the  particular  application.  For  the  case  of 
fractal  trees  the  fractal  dimension  Df  is  the  measure  of  the  space  that  a 

self-similar  structure  fills,  and  it  varies  with  the  branching  levels  that  the  structure 
consists  of. 

The  most  useful  terms  from  fractal  geometry  that  are  used  to  describe  a  planar 
fractal  tree  are  the  following: 

1.  The  number  of  branch  segments  N  formed  from  each  preceding  branch 
segment. 

2.  The  constant  similarity  ratio  "r"  that  relates  the  fractional  reduction 
in  segment  length  for  each  segment  to  previous  level,  this  factor  is  less 
than  1.0. 

In  this  case,  the  fractal  dimension  Df  is  given  by  the  formula  [Ref.  7]: 

1  -  Df  =  log( total  branch  length)  _  log(rN) 
log(average  branch  length)  log( 1 / r) 
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This  equation  is  accurate  for  symmetric  structures  or  asymmetric  structures  having 
a  normal  distribution  function  [Ref.  7].  More  details  about  fractal  geometry  are 
discussed  in  [Ref.  1]  and  [Ref.  7]. 

B.  FRACTAL  TREE  GENERATION 

In  the  planar  fractal  structures  that  are  used  for  calculating  the  radar  cross 
section,  the  number  of  branch  segments  N  formed  from  each  preceding  branch 
segment  is  2.  The  principal  properties  of  these  trees  are  that  the  branches  are  not 
overlapping,  and  the  angle  6  between  any  two  branches  is  the  same.  The 
nonoverlapping  between  the  branches  of  each  fractal  model  is  achieved  by  choosing 
the  branch  angle  $  for  a  given  reduction  factor  by  the  following  empirical  formula 
[Ref.  7] 


Minimum  0  -  32.34  x  r®^7  (Radians) 

This  angle  is  given  in  radians  and  r  is  the  desired  reduction  factor  which  is  a 
constant  for  the  structure.  As  the  reduction  factor  r  increases  the  tree  fills  more 
space. 

Five  types  of  planar  fractal  trees  are  considered  in  this  thesis.  Each  type 
corresponds  to  a  set  of  values  for  reduction  factor  r  and  minimum  branch  angle  0. 
Table  4.1  shows  these  sets  of  r  and  minimum  6,  win  re  0  is  given  in  degrees. 

Each  model  is  generated  by  selecting  the  desirable  reduction  factor  r  (r  <  1.0), 
the  corresponding  branch  angle  6  given  by  equation  4.1,  the  initial  length  from 
which  the  reduction  will  start,  the  value  of  N  (which  in  the  investigated  case  is  2), 
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and  the  number  of  the  desired  levels  for  branch  generation.  Each  level  contains  2r 
points  corresponding  to  the  end  points  of  the  reduced  length  branches. 


TABLE  4.1 

NUMERICAL  VALUES  OF  r  AND  MINIMUM  ANGLE  9. 


REDUCTION  FACTOR  frl  ANGLE  9 

0.53  29.940 

0.55  46.430 

0.60  95.890 

0.66  145.350 

0.71  180.000 


Figures  4. 1-4.5  show  the  planar  fractal  trees  that  correspond  to  each  set  of  r 
and  9  of  Table  4.1  respectively.  As  the  branch  angle  9  increases,  the  spreading  of  the 
branches  becomes  larger.  In  all  trees  the  physical  length  of  the  initial  branch  is 
chosen  as  two  centimeters.  These  models  were  generated  by  a  FORTRAN  program 
which  has  the  following  input  data: 

1.  The  initial  length.  This  the  length  of  the  first  dipole  whose  length  will 
be  reduced  by  the  constant  reduction  factor  (r). 

2.  The  value  of  N,  being  2  in  all  cases. 

3.  The  initial  point  that  the  reduction  starts.  This  poiut  is  the  end  point 
of  the  initial  wire. 
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4.  The  number  of  desirable  levels  for  branch  generation. 

5.  The  branch  angle  0. 

The  program  generates  two  output  files.  The  first  is  the  input  file  for  the  RCS 
program  containing  all  the  required  data  that  were  mentioned  in  Chapter  3.  The 
second  output  file  contains  the  coordinates  of  the  start  and  end  points  of  each 
generated  branch. 

For  the  models  characterized  by  reduction  factors  0.53,  0.55,  and  0.60,  the 
program  limits  the  structure  size  when  the  length  of  a  branch  becomes  smaller  than 
0.1  of  the  wavelength  of  the  incident  plane  wave.  For  the  models  characterized  by 
reduction  factor  0.66  and  0.71,  this  limit  is  taken  as  0.15.  The  reason  is  that  as  the 
reduction  factor  increases  the  number  of  branches  in  a  particular  branch  level 
increases.  The  size  of  the  tree  is  truncated  so  that  the  number  of  unknowns  is 
manageable. 

The  geometry  of  all  planar  fractal  trees  is  in  the  same  coordinate  system  that 
the  RCS  program  uses.  All  models  have  been  generated  in  the  y-z  plaue  of  a 
cartesian  coordinate  system. 
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Figure  4.3  Fractal  Tree  for  r  =  0.G0  and  0  =  95.89°. 
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Figure  4.5  Fractal  Tree  for  r  =  0.71  and  6  =  ISO0. 
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V.  NUMERICAL  RESULTS 


In  this  chapter,  computed  results  for  the  radar  cross  section  o  of  planar 
structures  are  presented.  In  order  to  check  the  program,  the  radar  cross  section 
results  of  a  single  planar  dipole  are  compared  with  the  results  given  in  [Ref.  8]. 
Details  of  this  comparison  are  presented  in  the  following  section. 

In  all  models  where  numerical  results  for  the  radar  cross  section  are  presented, 
the  following  factors  have  been  considered: 

1 .  E-plane  is  the  plane  corresponding  to  tpo  =  0°,  and  8q  varying  from  0°  to 

1800. 

2.  H-piane  is  the  plane  corresponding  to  6q  =  90°,  and  <po  varying  from  -90° 
to  +900. 

3.  The  resulting  RCS  a  is  normalized  to  the  square  of  the  wavelength  A  of  the 
incident  plane  wave  (<7/ A2). 

4.  The  number  of  segments,  that  each  dipole  is  subdivided,  is  taken  as  four 
per  wavelength. 

5.  The  physical  dimensions  and  orientations  of  the  dipoles  used  to  construct  a 
planar  structure  are  the  same  when  this  structure  is  investigated  at  different 
frequencies,  and  only  the  segmentation  is  different  depending  upon  the  wavelength  A 
of  the  incident  plane  wave. 

6.  Each  PWS  mode  has  length  equal  to  the  length  of  two  segments. 

7.  The  incident  electric  field  is  linearly  polarized  and  the  parameters  a  and  b 
are  taken  as  0  and  1  respectively  for  this  investigation.  This  corresponds  to  having 
only  a  horizontal  magnetic  field. 
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A.  SCATTERING  FROM  A  SINGLE  DIPOLE 


In  Figure  5.1  a  centered  loaded  vertical  planar  thin  dipole  of  length  L  and 

width  2w  is  oriented  along  the  z  axis.  This  dipole  is  excited  by  a  plane  wave  of 

frequency  30  GHz  traveling  in  the  x-y  plane  =  90°,  tpo  =  0°).  The  same 

situation  is  described  in  [Ref.  8,  pp.  510-515]  for  a  cylindrical  dipole  of  radius  a  and 

length  L  such  that  L/2a=74.2.  In  the  present  case  of  a  planar  dipole,  the  length  of 

the  planar  dipole  is  selected  such  that  L/2w  =  33  [Ref.  9].  Figure  5.2  shows,  on  a 

semi-log  scale,  the  results  computed  by  the  developed  program  of  the  monostatic 

o 

normalized  radar  cross  section  a/X  of  this  single  dipole  for  values  of  L/A  ranging 
from  0  to  1.4,  and  for  loads  Zl  =  0  and  Zl  =  oo.  The  case  of  Zl  =  oo  is  achieved  by 
setting  a  small  gap  (0.01  wavelength)  at  the  center  of  this  planar  dipole.  In  Figure 
5.2,  the  corresponding  values  of  a/X  from  [Ref.  8,  pp.  115]  are  also  shown. 

This  comparison  leads  to  an  accurate  check  of  the  correct  calculation  of  the 
radar  cross  section  by  the  developed  RCS  program.  The  small  discrepancy  between 
the  computed  values  and  those  given  by  [Ref.  8]  is  due  to  the  error  incurred  in 
reading  values  off  the  curves  presented  in  [Ref.  8,  pp.  115]. 
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Figure  5.1  Vertical  Cen 
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B.  SCATTERING  FROM  PLANAR  FRACTAL  TREES 

In  Chapter  4,  five  types  of  planar  fractal  trees  were  described.  Each  type  of 

these  trees  was  lying  in  the  y-z  plane.  The  developed  RCS  program  calculates  the 

broadside  radar  cross  section  a  of  each  tree  for  frequencies  15  GHz-75  GHz,  for  the 

monostatic  case.  The  initial  dipole  that  is  being  reduced  by  the  reduction  factor  r 

has  physical  length  2  cm  and  only  the  branch  lengths  are  changed  for  each  tree 

depending  upon  the  reduction  factor  r.  As  the  frequency  of  the  plane  wave  changes, 

the  electrical  length  of  this  dipole  as  well  as  the  dipoles  composing  the  tree  will  also 

change.  For  the  frequency  range  of  15  GHz-75  GHz,  the  electrical  length  of  the 

initial  length  will  vary  from  one  to  four  wavelengths. 

o 

Figures  5.3-5.6  show  the  variations  of  aj\  ,  in  dB,  in  the  E-plane  and  the 
H-plane  for  the  case  of  a  fractal  tree  characterized  by  reduction  factor  r  =  0.53  and 
branch  angle  $  =  29.94°.  The  frequency  of  the  incident  plane  wave  varied  from  15 
GHz  to  60  GHz  in  steps  of  15  GHz.  The  tree  is  composed  of  31  dipoles.  This  number 
is  small  as  the  reduction  factor  is  small  and  the  branch  lengths  are  reduced  to  very 
small  values  at  low  levels. 

The  a/ A  variations  in  the  E-Plane  are  in  a  range  of  10  dB  approximately.  In 
2 

the  H-Plane.  The  c/X  is  symmetric  about  the  90°  axis  and  is  smoother  at  lower 
frequencies.  As  the  frequency  increases,  the  maximum  value  of  a/ A  at  0q  =  90°  and 
po=0°,  increases  from  2.73  dB  at  15  GHz  to  19.23  dB  at  60  GHz.  Figure  5.7  shows 
the  variation  of  the  maximum  aj  A  in  terms  of  frequency  increments  for  the  same 
structure.  This  variation  is  very  small  between  30  and  45  GHz. 
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F  =  30  GHz 
E-Plane 


Figure  5.4  Normalized  RCS  at  F  -  30  GKz  of  a  Fractal  Tree 
with  r  =  0.53  and  6  -  29.94°. 
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Figure  5.5  Normalized  RCS  at  F  =  45  GHz  of  a  Fractal  Tree 
with  r  -  0.53  and  0  =  29.94°. 
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F  =  60  GHz 
E-Plane 


Figure  5.6  Normalized  RCS  at  F  =  60  GHz  of  a  Planar  Fractal  Tree 
with  r  =  0.53  and  0  =  29.940. 
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Figure  57  Variation  of  maximum  cjX  vs.  Frequency 
of  a  Flan  or  Fractal  Tree  with  r  =  0.53  and  9  =39.94°. 
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Figures  5.8-5.12  show  the  variations  of  a/X  ,  in  dB,  in  the  E  and  H  planes 
for  the  fractal  tree  characterized  by  a  reduction  factor  r  -  0.55  and  branch  angle  9 
=  46.43°.  In  this  case,  the  fractal  tree  is  composed  of  63  dipoles  and  spreads  more 
than  the  previous  one,  and  the  lengths  of  the  branches  are  bigger  than  the  previous 
ones  (due  to  a  larger  reduction  factor  of  0.55  instead  of  0.53).  This  results  in  more 

o 

variations  for  a/ X  in  both  E  and  H  planes. 

This  model  is  investigated  at  a  frequency  range  of  15—75  GHz.  At  all 
o 

frequencies,  the  a/X  varies  in  a  range  of  10-12  dB  approximately.  In  the  H-Plane 
2 

the  variation  of  a/X  is  symmetric  about  the  90°  axis.  The  maximum  value  of  radar 
cross  section  varies  from  2.78  db  at  15  GHz  to  22.63  dB  at  75  GHz.  Figure  5.13 

9 

shows  the  maximum  value  of  a/X  ,  corresponding  to  $o  =  90°  and  ipo  =  0°  in  terms 
of  the  frequency  of  the  incident  plane  wave,  varying  from  15  GHz  to  75  GHz.  This 
maximum  a/X  increases  linearly  as  the  frequency  increases  except  the  range  of 
45-60  GHz  where  the  variation  is  very  small. 
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Figure  5.S  Normalized  RCS  at  F  =  15  GHz  of  a  Planar  Fractal 
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Figure  5.9  Normalized  RCS  at  F  =  30  GHz  of  a  Planar  Fractal  Tree 
with  r  s=  0.55  and  6  ~  46.43°. 


Figure  5.10  Normalized  RCS  at  F  =  45  GHz  of  a  Planar  Fractal  Tree 
with  r  -  0.55  and  0  =  46.43°. 
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Figure  5. 11  Normalized  HCS  at  F  —  60  GHz  of  a  Flanar  Fractal  Tree 
with  l  -  0.55  and  0  -  46.43°. 
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Figure  5.12  Normalized  RCS  at  F  =  75  GHz  of  a  Planar  Fractal  Tree 
with  r  -  0.55  and  $  -  46.43°. 


Figures  5.14-5.18  show  the  variations  of  the  normalized  radar  cross  section 

2 

cr/A  of  a  planar  fractal  tree  characterized  by  reduction  factor  r  =  0.60  and  branch 
angle  6  =  95.89°.  This  model  is  composed  of  63  dipoles  and  has  more  spreading  than 
the  two  models  that  were  previously  investigated.  The  branches  have  larger  physical 
lengths  than  the  ones  in  other  two  models.  At  15  GHz,  the  variation  of  aj\  in  the 
E-Plane  is  very  low  and  similar  to  the  variation  in  the  other  two  models  at  the 
same  frequency.  In  the  H-Plane  this  variation  is  different  as  the  incident  plane 
wave  strikes  more  dipoles  from  -90°  to  90°. 

As  the  frequency  increases,  the  electrical  length  of  the  dipoles  is  also  increased, 

and  more  lobes  appear  at  30  GHz  and  higher  frequencies  in  both  the  E  and  the  H 

planes.  The  maximum  value  of  aj A  at  0o=  90°  and  tpo  =  0°  varies  from  -0.23  dB 

at  15  GHz  to  21.60  dB  at  75  GHz.  Figure  5.19  shows  the  variation  of  the  maximum 

radar  cross  section  in  terms  of  the  frequency  of  the  incident  plane  wave.  In  a  range 

o 

of  15-60  GHz  the  values  of  maximum  a/ A  follow  a  straight  line  approximately.  In 
the  range  of  60-75  GHz  the  variation  of  maximum  cr/A  is  small. 

The  planar  fractal  trees  characterized  by  r  =  0.66,  0  -  145.35°,  and  r  =  0.71,  0 
-  180°,  have  large  values  of  branching  angles.  The  large  values  of  reduction  factor  r 
generate  fractal  trees  whose  branches  have  large  physical  lengths  compared  with  the 
lengths  of  the  brandies  of  the  previously  investigated  fractal  trees.  The  result  is  that 
the  electrical  lengths  of  these  branches  are  large  also  at  the  range  of  15-75  GHz,  and 
a  large  number  of  modes  is  required  to  investigate  the  last  two  types  of  fractal  trees. 
It  was  found  that  the  developed  RCS  program  is  not  able  to  calculate  the  radar 
cross  section  of  fractal  trees  that  are  characterized  by  large  values  of  reduction 
factor  r  and  branching  angle  0  due  to  memory  restrictions  and  other  numerical 
problems. 
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2 

For  this  reason  the  numerical  results  for  a/X  are  not  presented  in  this  thesis  for 
these  two  cases. 

It  is  seen  by  comparing  Figures  5.7,  5.13,  and  5.19  that  there  is  a  frequency 
range  over  which  the  variation  of  maximum  RCS  is  rather  small.  Furthermore,  this 
range  of  frequencies  shifts  to  higher  frequencies  as  the  tree  structure  is  spread  from  r 
=  0.53  to  r  =  0.60.  The  maximum  RCS  of  each  of  these  trees  varies  in  a  range  of  3 
dB  approximately  at  the  same  frequency. 

Scattering  from  an  actual  tree  will  not  exactly  follow  all  these  patterns,  but  it 
is  felt  that  the  trends  would  generally  remain  the  same.  This  is  specially  true  for  the 
frequency  behavior. 
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Figure  5.15  Normalized  RCS  at  F  =  30  GHz  of  a  Planar  Fractal  Tree 
with  r  =  0.60  and  0  =  95.S90. 
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Figure  5.16  Normalised  RCS  at  F  -  75  GHZ  of  a  Planar  Fractal  Tie? 
with  r  =  0.60  and  0  -  95.60°. 


01 

o  cQ 

tC  LO 

o 

ii  ii 
© 


id  o  m 

CM  CM  ^ 


o 

oi 


in 

N 


o 


in 


o 

co 


in 


o 


(gp)  y  /  jd 

z 


Hgure  5.19  Variation  of  Maximum  RCS  vs.  Frequency 
of  a  Planar  Fractal  Tree  with  r  =  0.00  and  0  -  95.80°. 
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VI.  CONCLUSIONS  AND  RECOMMENDATIONS 


The  radar  cross  section  of  natural  trees  is  a  complicated  problem  of 
electromagnetic  scattering.  The  real  trees  are  made  up  of  long,  intersecting,  and 
lossy  objects.  The  geometry  of  these  objects  is  not  easy  to  set  up  in  a  three 
dimensional  space. 

in  this  thesis,  the  real  fractal  trees  were  approximated  by  planar  fractal  trees. 
The  radar  cross  section  of  these  trees  was  calculated  by  developing  a  Fortran 
program.  The  planar  fractal  trees  that  were  used  for  this  investigation  were 
symmetric  planar  structures  composed  of  perfectly  conducting  and  non-intersecting 
planar  dipoles.  The  geometry  of  these  structures  was  generated  using  fractal  theory. 

A.  CONCLUSIONS 

The  radar  cross  section  of  the  planar  fractal  trees  was  calculated  using  the 
moment  method  theory.  For  a  given  planar  structure  composed  of  planar  strips  and 
illuminated  by  a  plane  wave,  the  program  generates  the  geometry  of  the  P\V$ 
modes,  calculates  the  impedance  between  any  two  of  them,  and  fills  the  imjredance 
matrix  (Zj.  The  voltage  on  eacii  PWS  mode,  due  to  the  incident  electric  field,  is  also 
calculated,  and  the  voltage  column  vector  (Vj  is  generated. 

The  RCS  program  uses  the  moment  method  theory  to  calculate  the  current 
distribution  on  each  PWS  mode.  The  radiation  equations  of  electromagnetic  theory 
are  used  to  determine  the  scattered  electric  field  due  to  the  current  induced  on  each 
PW  S  mode  of  the  structure.  The  knowledge  of  both  the  given  incident 
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electric  and  calculated  scattered  electric  fields  leads  to  the  calculation  of  the  radar 

o 

cross  section  (<r/A  ),  for  the  given  angles  0o,<A)  of  the  incident  electric  field. 

The  calculated  radar  cross  section  of  a  single  centered  loaded  vertical  dipole 
was  compared  with  the  values  given  in  [Ref.  8,  pp,  115],  for  a  similar  dipole. 
Although  the  investigation  was  for  a  planar  dipole  and  in  [Ref.  8]  the  dipole  was 
cylindrical,  the  discrepancy  of  the  results  was  very  small. 

In  this  thesis,  five  models  of  planar  fractal  trees  were  investigated.  The 
geometry  of  these  models  was  generated  by  a  given  program  which  was  modified  to 
generate  the  input  data  for  the  developed  RCS  program. 

The  broadside  radar  cross  section,  for  the  monostatic  case,  was  calculated  for 
three  of  these  planar  fractal  trees,  for  a  frequency  range  of  15  GHz  to  75  GHz  in 
steps  of  15  GHz.  This  investigation  showed  that  the  radar  cross  section  varied  in  a 
range  of  10  dB  in  the  E-Plane.  In  the  H-Plane,  the  variation  of  ojX  was 
symmetric  about  the  90°  axis  and  smooth  at  low  frequencies  and  small  branching 
angles.  Fo:  higher  branching  angles  and  reduction  factors,  more  variations  of  oj\~ 
with  the  frequency  were  seen  in  both  the  E  and  II  planes. 

The  maximum  RCS  at  0$  =  90°,  tpo  =  0°  showed  a  small  variation  over  a 
frequency  range.  This  range  was  different  on  each  investigated  tree.  The  variation  of 
the  maximum  RCS  followed  a  straight  line  at  the  other  frequencies.  In  each  tree  and 
at  the  same  frequency,  the  maximum  RCS  varied  in  the  range  of  3  dll 
approximately. 

It  was  found  that  the  developed  RCS  program  was  not  able  to  calculate  the 
radar  cross  section  of  fractal  trees  that  were  characterized  by  large  values  of 
reduction  factor  and  branching  angle  due  to  memory  restrictions  and  other 
numerical  problems. 


65 


The  scattering  from  an  actual  tree  will  not  exactly  follow  the  patterns  that 
were  described  in  Chapter  5,  but  it  is  felt  that  the  trends  would  generally  remain 
the  same.  This  is  specially  true  for  the  frequency  behavior. 

B.  RECOMMENDATIONS 

One  recommendation  is  to  investigate  the  radar  cross  section  of  planar  fractal 
trees  characterized  by  values  other  than  those  used  in  this  thesis.  Especially,  a 
limited  set  of  reduction  factor  and  branching  angle  has  to  be  established  for  the 
same  physical  length  of  the  initial  dipoie. 

Another  recommendation  is  to  investigate  the  radar  cross  section  of  planar 
fractal  trees  characterized  by  the  same  values  of  reduction  factor  and  branching 
angles  as  those  that  were  used  in  this  thesis  but  with  different  physical  and 
electrical  dimensions  of  the  fractal  trees. 

In  this  thesis,  the  branches  of  the  fract  d  trees  were  assumed  to  be  perfectly 
conducting  planar  strips.  The  trees  were  considered  without  leaves.  The  radar  cross 
section  of  fractal  trees  composed  of  lossy  planar  strips  and  with  leaves  should  be 
investigated. 

Finally,  the  generation  of  a  three  dirneusional  fractal  tree  and  its  radar  cross 
section  may  be  investigated. 
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APPENDIX  A 
PROGRAM  RCS 


THIS  PROGRAM  CALCULATES  THE  RADAR  CROSS  SECTION  OF  A 
PLANAR  STRUCTURE  COMPOSED  OF  ARBITRARILY  ORIENTED 
PLANAR  THIN  DIPOLES.  THE  DIPOLES  ARE  NOT  INTERSECTING. 
THE  DIPOLES  ARE  PERFECT  CONDUCTORS. 


GEOMETRY  IN  Y-Z  PLANE 

INPUT  DATA: 

F  =  FREQUENCY  IN  GHZ 
NW  =  NUMBER  OF  DIPOLES 

A,B  =  COORDINATES  OF  INCIDENT  ELECTRIC  FIELD  ON  Y  AND  Z 
AXIS  RESPECTIVELY. 

LW  =  HALF  LENGTH  OF  EACH  DIPOLE  IN  CM 

WW  =  HALF  WIDTH  OF  EACH  DIPOLE  IN  CM 

SW,  TW  =  COORDINATES  OF  THE  CENTER  OF  EACH  DIPOLE 

ALONG  Y  AND  Z  AXIS  RESPECTIVELY 

PSIW  =  ANGLE  IN  DEGREES  BETWEEN  THE  Z  AXIS  (REFERENCE) 
AND  DIRECTION  OF  CURRENT  FLOW  ON  EACH  DIPOLE 
(POSITIVE  ANGLES  ARE  MEASURED  CLOCKWISE) 

NW  =  NUMBER  OF  SEGMENTS  THAT  EACH  DIPOLE  IS 
SUBDIVIDED 

THITAO,  PHIO  =  ANGLES  IN  DEGREES  OF  INCIDENT  ELECTRIC 
FIELD. 

OUTPUT  DATA: 


(' 

C 


NORMALIZED  RCS  (a/r)  IN  dB 

REAL  LS(  1 :230),S(  1:230), T(  1 :23Q),PSI(1 :230) 

REAL  LEN,PG.SU,TU,WID,A,B,MAG(  1:230) 

INTEGER  I,NMAX,INDX(230),NT,NP,L,M,N,K,NW 
COMPLEX  Z(i :230,1 :230),V(  1:230), CUR(4:230) 

COMMON/  RCI /  NMAX 
COMMON/  RC2/  LS,WI,S,T,PSI 
COMMON/  RC4/  CUR 

OPEN  (UNrr-l.EILE^'INI'UTr.EORM-TORMA'n'EI)') 
OPEN  ( UN IT-2,EI LE-'OU  I  PUT 1 ' ,FORM = ' UN  FORM  ATTED') 
OPEN  (UNIT=3.FILEa^,INPUT^FORM=,FORMA,^^5I),) 
OPEN(UNTT=4,FILE=,OUTPUT',FORM=T'ORMATTED’) 
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READ(1,*)  F,NW,A,B 
PI  =  4.*  ATANfl.) 

KO  =  PI*F/15. ' 

RAD  =  PI/180. 

CALL  ZMATR  (F,NW,Z) 

NP=230 

NT=NMAX 

CALL  CLUDCP  (Z,NT,NP,INDX) 

READ  (3,*,END=11)  THITA0,PHI0 
CALL  VOLT  F,THITAO,PHIO,A,B,V) 

CALL  CLUBSB  (Z,NT,NP,INDX,V) 

DO  54  K=1,NT 
CUR(K)=V(K) 

CONTINUE 

CALL  RADIAT  (F,THITAO,PHIO,NMAX,A,B,RCS) 

WRITE(4,*)  THITAO,PHIO,RCS 
PRINT*,  RCS 
GOTO  22 
CLOSER) 

STOP 

END 

SUBROUTINE  TO  COMPUTE  THE  MUTUAL  IMPEDANCE  BETWEEN 
THE  PWS  MODES  OF  N  ARBITRARILY  ORIENTED  DIPOLES. 

GEOMETRY  IN  THE  Y-Z  PLANE 

SUBROUTINE  ZMATR  (F,NW,Z) 

REAL  F,S(1 :230),T(  1 :230),LG,SG,TG, WI(  1 :230) ,PSI(  1 :230) 

REAL  PSIG, PSIW(1:70),LW(1:70),H,W,DH,DW, WIDTH, LENG 
REAL  H 1  ,H2,Wl ,W2,Sl ,S2,T1  ,T2,PSI  1  ,PS12,LS(  1 :230) 

REAL  W W(  1 :70),SW ( 1 :70) ,TW ( 1 :70),SS(  1 :20),TS(  1 :20) 

REAL  PG,WID,LEN,SU.TU,A,B 
INTEGER  TEMP, P, L,NS(1:70  , NMAX, I, N,NW 
INTEGER  FUN,G,J,GB,GR,K 
COMPLEX  Z12,Z(  1:230, 1:230) 

COMMON/  HC2/  LS,WLS,T,PS1 

COMMON/  RC1  /  NMAX 

COMMON/  JOHN/  LG,NG,SG,TG,PSIG,LENG 

COMMON/  ZPAR /  H,W,DW,I)H 

COMMON/  ZINC/  II 1 ,112, W1  ,W2,S1  ,S2,TI  ,T2,PS1 1,PS!2 

DO  101  I  =  I,  NW 

HEAl)(I,*)  LW(I),WW(1),SW(I),TW(I),PSIW(I),NS(I) 

CONTINUE 

CLOSE(l) 

NMAX -  0 
DO  75  1=1, NW 

NMAX  =  NMAX  +  (NS  (I)  -  1) 

CONTINUE 
P—  1 
GB~0 
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DO  50  1=1, NW 

LENG  =  LW(I) 

LG=2*LW(I)/FLOAT(NS(I)) 

WIDTH=WW(I) 

NG=NS(I) 

SG=SW(I) 

TG=TW(I) 

PSIG=PSIW(I) 

CALL  GEOM(SS,TS) 
GR=NS(I)+GB-1 
M  =  P 

DO  60  J=1,NS(I)-1 
PSI(M)=PSIG 
LS(M)=LG 
WI(M)=WIDTH 
S  M)=SS(J) 
T(M)=TS(J) 

PRINT*, S(M),T(M) 

M  =  M  +  1 

60  CONTINUE 

GB=GR 
P=GR+1 

50  CONTINUE 
TEMP=NS(1)-1 
L=1 

G=NMAX 
DO  80  M=1,G 

IF(M.LE.TEMP)GOTO  85 
L=L+1 

TEMP=TEMP+NS(L)-1 

85  CONTINUE 

K=TEMP 
DO  90  N=M,G 

IF(N.LE.K)TIIEN 

H=LS(M) 

WaWI(M) 

DW=0 

DH=ABS(M-N)*H 
CALL  ZSDIP(F,Z12) 
Z(M,N)=Z12 
PRINT*, Z(M,N) 
WRITE(2)  Z(M,N) 
ELSE 
HI=LS(M) 

H2=L8(N) 

\Vl=WI(M) 

\V2=\V1(N) 

S1=S(M) 

S2=S(N) 

T1=T(M) 
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T2=T 
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PSI1=PSI(M) 
PSI2=PSI  N) 

CALL  ZPSUR  (F,Z12) 
Z  M,N)=Z12 
PRINT*, Z(M,N) 
WRITE(2)  Z(M,N) 
ENDIF 
CONTINUE 
CONTINUE 
DO  24  M=2,G 

DO  26  N=1,M— 1 

Z(M,N)=Z(N,M) 
WRITE(2)  Z(M,N) 
CONTINUE 
CONTINUE 
RETURN 
END 


SUBROUTINE  TO  COMPUTE  THE  MUTUAL/SELF  IMPEDANCE 
BETWEEN  TWO  DIPOLES.  THE  DIPOLES  ARE  ASSUMED  TO  BE 
COPLANAR,  IDENTICAL  AND  PARALLEL.  THE  DIPOLE  TO 
DIPOLE  IMPEDANCE  IS  COMPUTED  AS  THE  SUM  OF  FOUR 
MONOPOLE  TO  MONOPOLE  IMPEDANCES. 

INPUT  PARAMETERS: 

F  ==  Frequency  of  operation  (GHz) 

H  =  Half  height  of  the  dipole  (cm) 

W  -  Half  width  of  the  dipole  (cm) 

DH  =  Longitudinal  distance  between  the  two  dipoles  (cm) 

DW  =  Transverse  distance  between  the  two  dipoles  (cm) 

SUBROUTINE  ZSDIP  (F,Z12) 

REAL  H,  W,  DW,  DH,  F 
COMPLEX  Z12,  ZT 
COMMON/  ZPAR/  H,W,DW,DH 
ZT  «  (Q.,0.) 

Z12  =  (0.,0.) 

CALL  ZSMONP  (F,  II,  W,  DW,  DH,  0, 1, 0, 1,  ZT) 

Z12  =  Z12  +  ZT 

CALL  ZSMONP  (F,  II,  W,  DW,  DH  +  H,  0, 1, 1, 0,  ZT) 

Z12  «  Z12  +  ZT 

CALL  ZSMONP  (F,  H,  W,  DW,  DH  -H,  1, 0,  0, 1,  ZT) 

Z12  a  Z12  +  ZT 

CALL  ZSMONP  (F,  H,  W,  DW,  DH,  1,  0, 1, 0,  ZT) 

ZI2  =  Z12  +  ZT 

RETURN 

END 
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SUBROUTINE  TO  COMPUTE  THE  SELF/MUTUAL  IMPEDANCE 
BETWEEN  TWO  IDENTICAL,  COPLANAR  MONOPOLES.  THE 
CURRENT  IS  ASSUMED  TO  BE  CONSTANT  IN  THE 
TRANSVERSE  DIRECTION. 

REF:  R.  JANASWAMY,  A  SIMPLIFIED  EXPRESSION  FOR  THE 
SELF/MUTUAL  IMPEDANCE  BETWEEN  TWO  COPLANAR 
AND  PARALLEL  MONOPOLES,  IEEE  T-AP,  AP-35, 

No.  10,  pp.  1174-1176,  October  1987. 

INPUT  PARAMETERS: 

F  =  FREQUENCY  IN  GHz 
H  =  LENGTH  OF  EACH  MONOPOLE  (cm) 

W  =  WIDTH  OF  EACH  MONOPOLE  (cm) 

D  =  CENTER  TO  CENTER  SPACING  BETWEEN  THE  TWO 
MONOPOLES  IN  THE  DIRECTION  TRANSVERSE  TO  THE 
CURRENT  FLOW  (cm) 

HH  =  CENTER  TO  CENTER  SPACING  BETWEEN  THE  TWO 
MONOPOLES  IN  THE  DIRECTION  OF  CURRENT  FLOW  (cm) 

111  =  TERMINAL  CURRENT  OF  END  1  OF  MONOPOLE  1. 

121  =  TERMINAL  CURRENT  OF  END  2  OF  MONOPOLE  1. 

112  =  TERMINAL  CURRENT  OF  END  1  OF  MONOPOLE  2. 

122  =  TERMINAL  CURRENT  OF  END  2  OF  MONOPOLE  2. 


NOTE:  Ill,  121, 112, 122  can  assume  values  only  0  or  1. 

ICODE  =  0,  IF  D  .LE.  4W 
ICODE  =  1,  OTHERWISE 

With  ICODE  =  0,  the  expression  provided  in  the  above  paper  is  used. 
With  ICODE  =  1,  a  modified  form  of  the  expression  provided  in  the 
above  paper  is  used.  (cf.  notes) 

OUTPUT  PARAMETERS: 


212  =»  COMPLEX  IMPEDANCE  BETWEEN  THE  TWO  SURFACE 
MONOPOLES. 


SUBROUTINE  ZSMONP  (F,  II,  W,  D,  HIE  Ill,  121, 112, 122,  Z12) 
REAL  F,  II,  W,  D,  HH,  A,  B,  PI,  K0,  V.  UB,  UBP,  UA,  UAP 
REAL  KW,  KH.  Kl),  RC.  FR,  FI,  II,  12. 13, 14,  SI,  Cl 
REAL  UABP,  2R,  21.  X,  AA  (1).  BB  l).  SQ.XV,  UAB.  TINY 
INTEGER  Ill,  121. 112, 122.  M,  N,  NX,  HI,  ICODE 
COMPLEX  212.  AC  (-1:1.  -1:1),  El.  E2,  21,  J,  CMN,  E3,  El,  FAC 
EXTERNAL  FR,  FI 


COMMON  / PARAM/  N.  V,  KD,  A.  B,  ICODE 
El  (X)aCI  (ABS  (X))  -  ,1  *  SI  (X) 

SQXV  (X)  -  SQRT  (X  *X  +  V"  V) 

TINY  =  LETS 
PI  =  4.  *  ATAN  (1.) 

NX  =  1 


1 1 


KO  =  PI  *  F  /  15. 

KW  =  KO  *  W 
KH  =  KO  *  H 
KD  =  KO  *  D 
ICODE  =  0 

IF  (KD  .GT.  4.  *  KW)  ICODE  =  1 
KI  =  10 

C 

C  Note  that  with  this  choice  of  KI  accurate  results  are  found  in  the 
C  region  0  <  D  <  4  W,  and  for  D  >  20  W.  In  between  these  two  regions, 

C  accurate  results  are  found  with  KI  =  20.  Hence  the  above  choice  is 

C  good  only  if  this  code  is  used  for  values  of  D  satisfying  the  above 
C  inequalities. 

C 

FAC  =  CMPLX  (1.,  0.) 

A  =  KD  -  2.  *  KW 
B  =  KD  +  2.  *  KW 
IF  (ICODE  .EQ.  1)  GO  TO  2 
AA  (1)  =  A 
BB  (1)  =  B 
GO  TO  3 

2  AA  (1)  =  -2.  *  KW 
BB  (1)  =  2.  *  KW 

3  II  =  (121  *  COS  (KH)  -  Ill)  *  112 
12=  Ill -121*  COS  (KH)  *122 

13  =  121 -111*  COS  KH  *112 

14  =  (121 -Ill*  COS  (KH))  *  122 
J  =  CMPLX  (0.,1.) 

Z1  =  CEXP  (J  *  KH) 

AC  (-1,  -1)  =  12  +  II  /  Z1 
AC  (-1,1)  =  12 +  11  *Z1 
AC  (1,  -1)  =  13  - 14  *  Z1 
AC  1,1)  =  13-14 /Z1 

AC  0,  -1)  =  -  (AC  (-1,-1)  *  Z1  +  AC  (1,  -1)  /  Zl) 

AC  (0, 1)  =  -  (AC  (1, 1)  *  Zl  +  AC  (-1. 1)  /  Zl) 

RC  =  15.  /  (2.  *  SIN  (KH)  *  KW)  **  2 
Z12  =  (0.,  0.) 

DO  1  M  =  -1, 1 
DO  1  N  =  -1,  1,  2 
V  =  K0  ■  (HU  +  M  *  II) 

IF  (ICODE  .EQ.  0)  GO  TO  4 
FAC  =  CEXP  (J  '  N  *  V) 

CMN  =  CMPLX  (0.,0.) 

GO  TO  5 

4  UA  “  SQXV  (A)  +  N  "  V 
UAP  =  UA  -  2.  *  N  *  V 
IB  =  SQXV  (U)  +  N  •  V 
UBP  =  UD  -  2.  *  N  *  V 
UAH  =  SQXV  (KD)  +  N  *  V 
UABP  =  UAB  -  2.  *  N  “  V 
El  =  El  (UU) 
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E2  =  (0.,0.) 

IF  (ABS  (A)  .GT.  TINY)  E2  =  El  (UA) 

E3  =  (0.,0.) 

IF  (ABS  (KD)  .GT.  TINY)  E3  =  El  (UAB) 

CMN  =  0.5  *  (B  *  B  *  El  +  A  *  A  *  E2  -  2.  *  KD  *  KD  *  E3  + 
k  CEXP  (-J  *  UB)  *  (1.  +  J  *  UBP)  +  CEXP  (-J  *  UA)  * 

k  (1.  +  J  *  UAP)  -  2.  *  CEXP  (-  J  *  UAB)  *  (1.  +  J  *  UABP)) 

5  CALL  HABER  (NX,  AA,  BB,  FR,  KI,  ZR) 

CALL  HABER  NX,  AA,  BB,  FI,  KI,  ZI) 

Z12  =  Z12  +  AC  (M,  N)  *  (CEXP  (J  *  N  *  V)  *  CMN  +  (ZR  +  J  *  ZI) 
k  *  FAC) 

1  CONTINUE 
Z12  =  Z12  *  RC 
RETURN 
END 
C 

REAL  FUNCTION  FR  (XI,  NX) 

INTEGER  N,  NX,  CODE 

REAL  X,  V,  Tl,  KD,  A,  B,  XI  (NX),  TKW,  T2,  Cl,  SI,  TINY 
COMPLEX  El,  J 

COMMON  /PARAM/  N,  V,  KD,  A.  B,  CODE 
El  (X)  =  CI  (ABS  (X))-J*SI  (X) 

TINY  =  l.E-6 
X  =  XI  (1) 

J  =  CMPLX  (0.,1.) 

IF  (CODE  .EQ.  1)  GO  TO  1 
Tl  =  SQRT  (X  *  X  +  V  *  V) 

FR  =  COS  (Tl) 

IF  (ABS  (V)  .GT.  TINY)  FR  =  FR  *  (1.  -  N  8  V  /  Ti) 

IF  X  .LE.  KD)  FR  =  FR  *  A 
IF  (X  .GT.  KD)  FR  =  -FR  *  B 
GO  TO  2 

1  TKW  =  0.5  *  (B  —  A) 

T2  =  SQRT  ((KI)  +  X)*(KD  +  X)  +  V*  V) 

FR  a  HEAL  (El  (T2  +  N  *  V)) 

FR  «  FR  -  (TKW  -  ABS  (X)) 

2  END 


REAL  FUNCTION  FI  (XI,  NX) 

INTEGER  N\  NX.  CODE 

REAL  X.  V,  Tl.  KD.  A,  B,  XI  (NX).  TKW,  T2.  Cl,  SI  TINY 
COM  PI  F  V  FI  1 

COMMON  /PARAM/  N.  V.  KI),  A.  B,  CODE 
El  (X)~d  (ABS  (X))  -  J  *  SI  (X) 

X  =  XI  (!) 

TINY  a  l.K-6 
J  =  CMPLX  {()..  1.) 

T!  =  SQRT  (X  *  X  -J-  V  *  V) 

IF  (CODE  .EQ.  I)  GO  TO  1 
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FI  =  -SIN  (Tl) 

IF  (ABS  (V)  .GT.  TINY)  FI  =  FI  *  (1.  -  N  *  V  /  Tl) 

IF  (X  .LE.  KD)  FI  =  FI  *  A 
IF  (X  .GT.  KD)  FI  *  -FI  *  B 
GO  TO  2 

TKW  =  0.5  *  (B  -  A) 

T2  =  SQRT  ((KD  +  X)  *  (KD  +  X)  +  V  *  V) 

FI  =  AIMAG  (El  (T2  +  N  *  V) 

FI  =  FI  *  (TKW  -  ABS  (X)) 

2  END 

SUBROUTINE  TO  COMPUTE  THE  COORDINATES  OF  THE  PWS 
MODES  OF  EACH  DIPOLE  IN  THE  Y-Z  PLANE 


225 


C 

c 


SUBROUTINE  GEOM  (S,T) 

REAL  THI,H1,H2,Y,Z,SG,TG,LG,PSIG 
REAL  P,Q,S(1:20),T(1:20),LENG 
INTEGER  M,NG,K 

COMMON/  JOHN /  LG,NG,SG,TG,PSIG,LENG 
COSD(X)  =  COS(X*RAD) 

SIND(X)  =  SIN(X*RAD) 

PI=4.*ATAN(1.) 

RAD=PI/180. 

THI=90.-PSIG 

Hl=COSD(THI) 

H2=SIND(THI) 

Y=SG-{LENG*H1) 

Z=TG-(LENG*H2) 

P=LG*Hl 
Q=LG*H2 
S(1)=*Y+P 
T(1)=Z+Q 
DO  225  K=2,NG-I 

S(K)=S(I)+(K-1)*P 

T(K)-T(1)+(K-1)*Q 

CONTINUE 

RETURN 

END 


SUBROUTINE  ZPSUR  (F.Z12) 

REAL  HU  Wl,  H2,  W2,  PSI,  F.  YSTAR.  ZSTAR.  K0,  A I  (2),  HI  (2) 
REAL  C  (3),  D  (3),  RT,  ZT,  COT,  CSEC,  KOS.  X.  PI,  PSU,  PSI2 
REAL  SI,  Tl.  S2,  T2,  IUNTG,  IMINTG,  RESULT,  SIND,  COSD.Sl.C1 


INTEGER  NX.  KI 
COMPLEX  Z12,  J 
EXTERNAL  HINTG.  IMINTG 


COMMON/  ZINC/  m.H2.\VT.VV2,Sl.S2.Tl,T2,PSU.PS12 
COMMON  /PARAM2/  C.  I).  RT,  ZT,  CSEC.  COT,  KOS,  J,  KQ 
SIND  (X)  =  SIN  (X  *  PI  /  ISO  ) 
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COSD  (X)  =  COS  (X  *  PI  /  180.) 

DATA  AI,  BI  /2  *  -1..  2  *  1./ 

PI  =  4.  *  ATAN  (1.) 

PSI  =  PSI2  —  PSI1 

IF  (ABS  (PSI).LF,0.4)  PSI=  SIGN(1.,PSI)*0.4 
J  =  CMPLX  (0.,1.) 

NX  =  2 
Ki  —  3 

KO  =  PI  *  F  /  15. 

C(l)  =  l./SIN  (KO  *  HI) 

C  3  =  C  (1) 

C  2)  =  -2.  *  COS  (KO  *  HI)  *  C  (1) 

D  1)  =  1.  /  SIN  (KO  *  H2) 

D  3  =D(1) 

D  (2)  =  -2.  *  COS  (KO  *  H2)  *  D  (1) 

CSEC  =  1.  /  SIND  (PSI) 

KOS  =  COSD  (PSI) 

COT  =  KOS  *  CSEC 

YSTAR  =  (S2-S1 )  *  COSD  (PSI1)  -  (T2-T1)  *  SIND  (PSI1) 

ZSTAR  =  (S2-S1)  *  SIND  (PSI1)  +  (T2-T1)  *  COSD  PSIl  j 

RT  =  YSTAR  *  CSEC  -  H2 

ZT  =  YSTAR  *  COT  -  HI  -  ZSTAR 

CALL  HABER  (NX,  AI,  BI,  RINTG,  KI,  RESULT) 

Z12  =  RESULT 

CALL  HABER  (NX,  AI,  BI,  IMINTG,  KI,  RESULT) 

Z12  =  Z12  +  J  *  RESULT 
Z12  =  -3.75  *  Z12 
RETURN 
END 

* 

REAL  FUNCTION  SI  (XI) 

REAL  XI 

REAL  AF  (4).  BF  (4),  AG  (4),  BG  (4),  X,  X2,  X4,  X6,  XS,  FX,  GX, 
&  PI.  SON- 

DATA  AF  /  38.027264.  265.187033,  335.67732,  38.102495  / 

DATA  BF  /  40.021433,  322.624911,  570.23628,  157.105423  / 

DATA  AG  /  42.242855,  302.757865,  352.018498.  21.821899  / 

DATA  BG  /  48.196927,  482.485984,  1114.978885,  449.690326  / 

SI  a  0. 

IF  (XI  .EQ.  0.)  RETURN 
SON  =  +1. 

IF  (XI  XT.  0.)  SON  «  -1. 

X  a  ABS  (XI) 

X2  =  X  "  X 

IF  (X  .GE.  20.)  THEN 

FX  «  1.  /  X  •  U.-2.  /  X2) 

GX  a  i./X2*  (l.-6./X2) 

GO  TO  1 
END  IF 
X  J  =  X2  *  X2 
X6  =  X4  #  X2 
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X8  =  X4  *  X4 
IF  (X  .LT.  1.)  THEN 

oJ  =  Xl1;,"?2  /  18-  +  X4  /  600.  - X6  /  35280.  +  X8  /  326592 \ ; 

SI  =  SI  ¥  SGN 

ELSE 

FX  =  (1.  +  AF  (1)  /  X2  +  AF  (2)  /  X4  +  AF  (3)  /  X6  +  AF  (4)  / 

X8) 

Fx  =  nxr(lw *Jin+  BF  (1)  /  X2  +  BF  (2)  /  X4  +  BF  (3)  /  X6  + 

GX  =  ^8)+  AG  (1)  1  X2  +  AG  1  X4  +  AG  (3)  1  X6  +  AG  1 

GX  =  SS/Kvoi1*  +  BG  (1)  /  X2  +  BG  (2)  /  X4  +  BG  (3)  /  X6  + 

PI  =  4.  *  At/n  (1.) 

X  =  X  -  AINT  (X  /  (2.  *  PI))  *  2.  *  PI 

SI  =  SGN  *  (PI  /  2.  -  FX  *  COS  (X)  -  GX  *  SIN  (X)) 

END  IF  J 

END 


REAL  FUNCTION  Cl  (X) 

a  tL  (4)’  BF  (4)’  BG  (4)>  x>  X2> X4’  X6>  X8>  Fx>  GX 

REAL  PI 

DATA  AF  /  38.027264,  265.187033,  335.67732,  38.102495  / 

DATA  BF  /  40.021433,  322.624911,  570.23628,  157.105423  / 

DATA  AG  /  42.242855,  302.757865,  352.018498.  21.821899  / 

DATA  BG  /  48.196927,  482.485984,  1114.978885,  449.590326  / 

IF  (X  .LE.  0.)THEN 

PRINT  *,  'Invalid  argument  for  Cl  (x)\  '  x  =  ',  x 

RETURN 

ELSE 

X2  =  X  *  X 
Xi  =  X2  *  X2 


IF  (X  .GE.  20.)  THEN 
FX  ss  1.  /  X  *  (1.  -  2.  /  X2) 

GX  =  1.  /  X2  *  (I.  -  6.  /  X2) 

GO  TO  I 
END  IF 
X6  =  XI  *  X2 
X8  =  X4  *  X4 
IF  (X  .LT.  I.)  THEN 

Cl  =  0.57721566  4-  A  LOG  (X)  -  X2  *  (0.25  -  X2  /  90.  +  Xi  /  4320. 
-  XG  /  322560.) 

ELSE 

EX  ^  (U  AK  (1)  /  X2  +  AF  (2)  /  X4  +  AF  (3)  /  XG  +  AF  (4)  / 
Xb) 

Kx = BF  0 }  1  x2 + uf  (2)  - x4 + bf  (3)  / x6  * 
GX  «  (t  -f  AG  (I)/  X2  -f  AG  (2)  /  X4  +  AG  (3)  /  X6  +  AG  (4)  / 


GX  =  GX  /  (X2  *  (1.  +  BG  (1)  /  X2  +  BG  (2)  /  X4  +  BG  (3)  /  X6  * 
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&  BG(4)  /  X8)) 

1  PI  =  4.  *  ATAN  (1.) 

X  =  X  -  AINT  (X  /  (2.  *  PI))  *  2.  *  PI 
Cl  =  FX  *  SIN  (X)  -  GX  *  COS  (X) 

END  IF 
END  IF 
END 
C 

REAL  FUNCTION  RINTG  (X,  NX) 

COMPLEX  J,  El,  TERM1,  TERM2,  INTGR 

INTEGER  NX,  M,  N,  P,  Q 

REAL  *  4  PSI1,  PSI2,  SI,  Tl,  S2,  T2 

REAL  *  4  C  (3),  D  (3),  R,  Z,  RT,  ZT,  U,  V,  CSEC,  COT,  KOS,  Y 
REAL  *  4  X  (NX),  HI,  H2,  Wl,  W2,  KO,  PZQR,  SI,  C A 
&  ,TT,  ZM,  RN ,  RMN,  TESC,  TESD,  PI,  PZQRP 

COMMON/  ZINC/  H1,II2,W1,W2,S1,S2,T1,T2,PSI1,PSI2 
COMMON  /PARAM2/  C,  D,  RT,  ZT,  CSEC,  COT,  KOS,  J,  KO 
El  (Y)  =  Cl  (ABS(Y))  -  J  *  SI  (Y) 

U  =  X(1) 

V  =  X  (2) 

INTGR  =  (0.,0.) 

TESC  =  ABS  (C  (2)  /  C  (1)) 

TESD  =  ABS  (D  (2)  /  D  (1)) 

TT  =  ALOG  (ABS  ((l.+KOS)/(l -KOS))) 

PI  =  4.  *  ATAN  (1.) 

DO  3M  =  1,3 

IF  (M  .EQ.  2  .AND.  TESC  .LE.  l.E-6)  GO  TO  3 
Z  =  (M-l)  *  HI 

ZM  =  -U  *  WI  *  COT  4-  V  *  W2  *  CSEC  +  ZT  +  Z 
TERM2  =  (0.,0.) 

DO  2  N  -  1,  3 

IF  (N  .EQ.  2  .AND.  TESD  .LE.  l.E-€)  GO  TO  2 
R=(N-1)*H2 

RN  =  -U  *  Wl  *  CSEC  +  V  *  W2  *  COT  +  RT  +  R 
IF  (ABS  (KO  *  RN)  .LE.  0.5E-1)  THEN 
PZQRP  =  KO  *  ABS  (ZM)  -  AINT  (KO  *  ABS  (ZM)  /  (2.  *  PI))  * 
t  2.  *  PI 

TERM1  =  CEXP  (-J  *  PZQRP)  *  TT 
GO  TO  4 
END  IF 

IF  (ABS  (KO  *  ZM)  .LE.  0.5E-1)  THEN 
PZQRP  -  KO  *  ABS  (RN)  -  AINT  (KO  *  ABS  (RN)  /  (2.  *  PI))  * 
k  2.  *  PI 

TERM1  =  CEXP  (-J  *  PZQRP)  *  TT 
GO  TO  4 
END  IF 

RMN  =  SQRT  ''RN  *  RN  +  ZM  *  ZM  -  2.  *  RN  *  ZM  *  KOS) 
TERM1  =  (0.,0.) 

DO  1  P  =  -1,  1,  2 
DO  1  Q  =  -l,  1,  2 
PZQR  =  P  *  ZM  +  Q  *  RN 
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PZQRP  =  PZQR  *  KO  -  AINT  (PZQR  *  KO  /  (2.  *  PI))  *  2.  *  PI 
TERM1  =  TERM1  +  P  *  Q  *  CEXP  (J  *  PZQRP)  * 
k  El  (KO  *  (RMN  +  PZQR))  '  ' 

1  CONTINUE 

4  TERM2  =  TERM2  +  D  (N)  *  TERM1 

2  CONTINUE 

INTGR  =  INTGR  +  C  (M)  *  TERM2 

3  CONTINUE 

RINTG  =  REAL  (INTGR) 

END 

REAL  FUNCTION  IMINTG  (X,  NX) 

COMPLEX  J,  El,  TERM1,  TERM2, 'INTGR 

INTEGER  NX,  M,  N,  P,  Q 

REAL  *  4  PSI1,  PSI2,  SI,  Tl,  S2,  T2 

REAL  *  4  C  (3),  D  (3),  R,  Z,  RT,  ZT,  U,  V,  CSEC,  COT,  KOS,  Y 
REAL  *  4  X  (NX),  HI,  H2,  Wl,  W2,  KO,  PZQR,  SI,  Cl 
k  ,TT,  ZM,  RN,  RMN,  TESC,  TESD,  PI.  PZQRP 

COMMON/  ZINC/  H1,H2,W1,W2,S1,S2,T1,T2,PSI1,PSI2 
COMMON  /PARAM2/  C,  D,  RT,  ZT,  CSEC,  COT,  KOS,  J,  KO 
El  (Y)  =  Cl  (ABS(Y))  -  J  *  SI  (Y) 

U  =  X(1) 

V  =  X  (2) 

INTGR  =  (0..0.) 

TESC  =  ABS  (C  (2)  /  C  (1)) 

TESD  =  ABS  (D  (2)  J  D  (1)) 

TT  =  ALOG  (ABS  ((l.+KOS)/(l.-KOS))) 

PI  =  4.  *  ATAN  (1.) 

DO  3  M  =  1,  3 

IF  (M  .EQ.  2  .AND.  TESC  .LE.  l.E-6)  GO  TO  3 
Z  =  (M-l)  *  HI 

ZM  =  -U  *  Wl  *  COT  +  V  *  W2  *  CSEC  +  ZT  +  Z 
TERM2  =  (0.,0.) 

DO  2  N  =  1,  3 

IF  (N  .EQ.  2  .AND.  TESD  .LE.  l.E-6)  GO  TO  2 
R  =  (N-l)  *  H2 

RN  =  -U  *  Wl  *  CSEC  +  V  *  W2  *  COT  +  RT  +  R 
IF  (ABS  (KO  *  RN)  .LE.  0.5E-1)  THEN 
PZQRP  =  KO  *  ABS  (ZM)  -  AINT  (KO  *  ABS  (ZM)  /  (2.  *  PI))  * 
k  2.  *  PI 

TERM1  =  CEXP  (-J  *  PZQRP)  *  TT 
GO  TO  4 
END  IF 

IF  (ABS  (KO  *  ZM)  .LE.  .5E-1)  THEN 
PZQRP  =  KO  *  ABS  (RN)  -  AINT  (KO  *  ABS  (RN)  /  (2.  *  PI))  * 
k  2.  *  PI 

TERM1  =  CEXP  (-J  *  PZQRP)  *  TT 
GO  TO  4 
END  IF 

RMN  =  SQRT  (RN  *  RN  +  ZM  *  ZM  -  2.  *  RN  *  ZM  *  KOS) 


78 


TERM1  =  (0.,0.) 

DO  1  P  =  -l,  1,2 
DO  1  Q  =  -l,  1,  2 
PZQR  =  P  *  ZM  +  Q  *  RN 

PZQRP  =  KO  *  PZQR  -  AINT  (KO  *  PZQR  /  (2.  *  PI))  *  2.  *  PI 
TERM1  =  TERM1  +  P  *  Q  *  CEXP  (J  *  PZQRP)  * 
k  El  (KO  *  (RMN  +  PZQR)) 

1  CONTINUE 

4  TERM2  =  TERM2  4-  D  (N)  *  TERM1 

2  CONTINUE 

INTGR  =  INTGR  +  C  (M)  *  TERM2 

3  CONTINUE 

IMINTG  =  AIMAG  (INTGR) 

END 

C 

C  Subroutine  to  compute  a  sequence  of  estimates  EST1  (K)  and 
C  EST2  (K),  1  .LE.  K1  .LE.  K  .LE.  K2  for  the  N-dimensional  integral 
C  B1  BN 

C  Int ...  Int  FUN  (xl,  x2, ...  xN)  dxl  dx2...  dx3 

C  A1  AN 

C  by  Haber's  method. 

C  Ref:  P.  J.  Davis  and  P.  Rabinowitz,  Methods  of  Numerical 
C  Integratation,  Academic  Press,  1984. 

C 

C  For  each  estimate  EST1  (K),  two  additional  quantities  ERR1(K) 

C  and  DEVI  (K)  are  computed.  If  the  values  of  DEVI  (K)  do  not 

C  vary  by  more  than  10%  between  consecutive  values  of  K,  then 

C  ERR1  (K)  can  be  taken  as  a  reliable  bound  on  the  difference 

C  between  EST1  and  the  integral.  A  similar  situation  holds  for 

C  EST2,  DEV2,  and  EHR2  (K).  The  total  number  of  functional 
C  evaluations  is  4  *  (K1  **  N  +  (Kl+1)  **  N  +  ...  +  K2  **  N)  and  I<2 

C  should  be  chosen  so  as  to  make  this  may  be  halfed  by  eliminating  the 

C  computation  of  the  EST2  (K).  In  other  situations,  these  values 

C  are  much  better  than  the  EST1  (K).  A  program  FUNCTION  FUN  (X,  N) 

C  must  be  supplied  by  the  user  with  X  declared  by  the  statement 

C  DIMENSION  X  (N).  FUN  must  be  declared  EXTERNAL  in  the  calling 

C  program.  If  N  <  1  or  N  >  10  or  K1  <  1  or  K2  <  Kl,  the  program 

C  terminates  with  IND  =  0.  Otherwise  IND  =  1. 

C 

C  Modified  by  R.  Janaswamy  so  that  the  output  is  average  of  EST1 
C  EST2.  Also,  Kl  =  K2  =  K. 

C 

SUBROUTINE  HABER  (N,  LL,  UL,  FUN,  K,  RESULT) 

INTEGER  N,  IND,  KEY,  I,  K,  J 
DOUBLE  PRECISION  AL  (10),  BE  (10),  GA  (10),  B,  G 
REAL  FUN,  Yl,  Y2,  Y3,  Y4,  EST1,  EST2,  ERR1,  ERR2, 
k  DEVI,  DEV2,  RESULT 

REAL  *  8  SI,  S2,  Dl,  D2 

REAL  LL  (N),  UL  (N),  DEX  (10),  PI  (10),  P2  (10),  P3  (10),  P4  (10), 
k  Q1  (10),  Q2  (10),  Q3  (10),  Q4  (10),  RAN  (10),  AKN,  AK,  T,  JAC 
REAL  AK1,  BK 
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EXTERNAL  FUN 

DATA  AL/.4142135623730950,  .7320508075688773,  .2360679774997897, 
.6457513110645906, .3166247903553998, .6055512754639893, 
.1231056256176605, .3589989435406736, .7958315233127195, 
.3851648071345040/ 

IND  =  0 

IF  (N  .LT.  1  .OR.  N  .GT.  10)  RETURN 
IND  =  1 
JAC  =  1. 

DO  1 1  =  1,  N 
BE  (I)  =  AL  (I) 

GA  (I)  =  AL  (I) 

RAN  (I)  =  UL  (I)  -  LL  (I) 

JAC  =  JAC  *  RAN  (I) 

1  DEX  (I)  =  0. 

AK  =  FLOAT (K) 

KEY  =  0 
AK1  =  AK  —  1.1 

51  =  0. 

D1  =  0. 

52  =  0. 

D2  =  0. 

AKN  =  AK  **  N 
T  =  SQRT  (AKN)  *  AK 
BK  =  1.  /  AK 

5  KEY  =  KEY  +  1 

IF  (KEY  .EQ.  1)  GO  TO  6 
KEY  =  KEY  -  1 
J  =  1 

4  IF  (DEX  (J)  .GT.  AK1)  GO  TO  8 
DEX  (J)  =  DEX  (J)  +  1. 

GO  TO  6 

8  DEX  (J)  =  0. 

J  =■•  J  +  1 

IF  (J  .LE.  N)  GO  TO  4 
GO  TO  3 

6  DO  7  1  =  1 ,  N 

B  =  BE  (I)  +  AL  (I) 

IF  (B  .GT.  1.)  B  =  B-l. 

G  -  GA  (I)  +  B 

IF  (G  .GT.  1.)  G  -  G  - 1. 

BE  (I)  =  B  +  AL  (I) 

IF  (BE  (I)  .GT.  1.)  BE  (I)  =  BE  (I)  -  1. 

GA  (I)  =  BE  (I)  +  G 

IF  (GA  (I)  .GT.  1.)  GA  (I)  =  GA  (I)  - 1. 

PI  (I)  =  (DEX  (I)  +  G)  *  BK 
Q1  (I)  =  LL  (I)  +  RAN  (I)  *  PI  (I) 

P2  (I)  =  (DEX  (I)  +  1.  -  G)  *  BK 
Q2  (I)  =  LL  (I)  +  RAN  (I)  4  P2  (I) 

P3  (I)  =  (DEX  (I)  +  GA  (I))  *  BK 
Q3  (I)  =  LL  (I)  +  RAN  (I)  *  P3  (I) 
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P4(I)  =  (DEX(I)  +  1.-GA(1))*BK 
7  Q4  (I)  =  LL  (I)  +  RAN  (I)  *  P4  (I) 

Y1  =  FUN  (Ql,  N) 

Y2  =  FUN  (Q2.  N) 

Y3  =  FUN  Q3,  N 
Y4  =  FUN  (Q 4.  N) 

51  =  SI  +  Y1  +  Y2 

D1  =  D1  +  (Y1  -  Y2)  **  2 

52  =  S2  +  Y3  +  Y4 

D2  D2  +  (Y1  +  Y3  -  Y2  -  Y4)  **  2 
GO  TO  5 

3  EST1  =  0.5  *  SI  /  AKN 

ERR1  =  1.5  *  DSQRT  (Dl)  /  AKN 
DEVI  =  ERR1  *  T 
EST2  =  0.25  *  (SI  +  S2)  /  AKN 
ERR2  =  0.75  *  DSQRT  (D2)  /  AKN 
2  DEV2  =  ERR2  *  T  *  AK 

RESULT  =  0.5  *  (EST?  +  EST2)  *  ABS  MAC) 

RETURN 
END 

SUBROUTINE  VOLT  WHICH  CALCULATES  THE  VOLTAGE  MATRIX 
VM 

GEOMETRY  IN  THE  Y-Z  PLANE 


SUBROUTINE  VOLT(F,THITAO,PHIO,  A,  B,V) 

REAL  TSI,TCO,PS,PC,KY,KZ,A,B,KO,F,KX 

REAL  MAR,GTI,PRS,PRC,MOD,KZHTA 

REAL  S(1:230),T(1:230),PSI(1:230),LS(1:230),WI():230),TH1TA0,PHI0 

INTEGER  NMAX 

COMPLEX  V(1:230),J,STR,VOM 

COMMON/  BRAVO/  KX,KY,KZ 

COMMON/  RC2/  LS,WI,S,T,PSI 

COMMON/  RC1/  NMAX 

COSD(X)=COS(X*RAD) 

SIND(X)=SIN(X*RAD) 

PI=4.*ATAN(1.) 

RAD=PI/180. 

K0=PI*F/15. 

J=CMPLX(0.,1.) 

TSI=SIND(THITAO) 

TCO=COSD(THITAO) 

PS=SIND(PHI0) 

PC=COSD(PHIO) 

KY=K0*TSI*PS 
KX=K0*TSI*PC 
KZ=KO*TCO 
DO  234  K=1,NMAX 

PRS=SIND(PSI(K)) 
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PRC=COSD(PSI(K)) 

MAR=KY*S  K)+KZ*T(K) 

STR=CEXP  (-J*  MAR) 

VOM=(STR/SIN(LS(K)*KO))*(A*PRS+B*PRC) 
KZHTA=KY*PRS+KZ*PRC 
MOD=COS(KO*LS(K))-COS(KZHTA*LS(K)) 
GTI=KZHTA**2— K0**2 
IF(ABS(KZHTA-KO) /KO.LT.  1  .E-2.0R. 
k  .ABS(KZHTA+KO)/KO.LT.l.E— 2)THEN 

V(K)=VOM*LS(K)*SIN(KO*LS(K)) 

ELSE 

V(K)=(2*KO/GTI)*VOM*MOD 

ENDIF 

234  CONTINUE 
RETURN 
END 
C 

SUBROUTINE  CLUDCP  (A,  N,  NP,  INDX) 

INTEGER  NP,  N,  INDX  (NP).,  I,  J,  K,  P,  NMAX 
PARAMETER  (NMAX  =  230) 

COMPLEX  A  (NP,  NP),  TEMP,  ETA.  W  (NMAX) 
REAL  AAMAX,  DUM 
DO  1  K  =  1,  N— 1 
AAMAX  =  CABS  (A  (I\,  I<)) 

P  =  K 

DO  2  I  =  K+l,  N 
DUM  =  CABS  (A  (I,  K)) 

IF  (DUM  .GT.  AAMAX)  THEN 
AAMAX  =  DUM 
P  =  I 
END  IF 

2  CONTINUE 
INDX  (K)  =  P 
DO  3  J  =  1,  N 
TEMP  =  A  (K,  J) 

A  (K,  J)  =  A  (P,  J) 

A  (P,  J)  =  TEMP 

3  CONTINUE 

DO  4  J  =  K+l,  N 
W  (J)  =  A  (K,  J) 

4  CONTINUE 

DO  5  I  =  K+l,  N 

ETA  =  A  (I,  K)  /  A  (K,  K) 

A  (I,  K)  =  ETA 

DO  6  J  =  K+l,  N 

A  (I,  J)  -  A  (I,  J)  -  ETA  *  W  (J) 

6  CONTINUE 

5  CONTINUE 
1  CONTINUE 

RETURN 

END 
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SUBROUTINE  CLUBSB  (A,  N,  NP,  INDX,  X) 

INTEGER  N,  NP,  I,  J,  INDX  (NP),  NMAX 

PARAMETER  (NMAX  =  230 

COMPLEX  A  (NP,  NP),  X  (NP),  TEMP,  Y  (NMAX) 

DO  PERMUTATIONS  ON  THE  EXCITATION  VECTOR  USING  THE 
INFORMATION  ON  THE  ROW  OPERATIONS  DONE  IN  CLUDCP. 


DO  1  K  =  1,  N-l 
TEMP  =  X  (K) 

X  (K)  =  X  (INDX  (K)) 
X  INDX  (K))  =  TEMP 
CONTINUE 


FORWARD  ELIMINATION 


DO  2  I  =  1,  N 
Y  (I)  =  X  (I) 
DO  2  J  =  1,  1-1 


Y  (I)  =  Y  (I)  -  A  (I,  J)*Y(J) 
CONTINUE 


BACK  SUBSTITUTION 


DO  31  =  N,  1,-1 
X  (I)  -  Y  (I) 

DO  4  J  =  1+1,  N 

X  (I)  =  X  (I)  -  A  (I,  J)  *  X  (J) 

CONTINUE 

X  (I)  =  X  (I)  /  A  (I,  I) 

CONTINUE 

RETURN 

END 


SUBROUTINE  RADIAT  TO  COMPUTE  THE  RADAR  CROSS  SECTION 
OF  A  PLANAR  STRUCTURE  COMPOSED  OF  ARBITRARILY 
ORIENTED  PLANAR  DIPOLES. 

SUBROUTINE  RADIAT  (F,THITA0,PHI0,NMAX  \,B,RCS) 

REAL  RC, THITAO, PHIO, LS(1:230), W1(1:230),S(  1:230), T(l:230), PI 
REAL  CR,F,K0,PI,UN,DM,EM,rAd,P,R,A,B,THITA,PHI,RCS,PSI(1:23O) 
REAL  KX,KY,KZ,KRON,KG 
INTEGER  G, NMAX 

COMPLEX  J,NTHI0,NPHI0, TEMPI, TEMP2,CUR(1:230) 

COMMON/  BRAVO/  KX.KY.KZ 
COMMON/  RC2/  LS,WI,S,T,PSI 
COMMON/  RC4/  CUR 
COSD(X)  =  COS(X*RAD) 

SIND(X)  =  SIN(X*RAD) 
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PI  =  4.  *  ATAN(1.) 

RAD  =  PI/180. 

KO  =  PI*F/15. 

UN=120.*PI 
THITA  =180.—  THITAO 
PHI  =  180.  +  PHIO 
J=CMPLX(0.,1.) 

TEMP1=(0.,0.) 

TEMP2=(0.,0.) 

NTHI0=(0.,0.) 

NPHI0=(0.,0.) 

DO  58  G— 1,NMAX 
P1=PSI(G) 

DM=KO*(S(G)*SIND(THITA)*SIND(PHI)+T(G)*COSD(THITA)) 
EM=K0*(SIND(P1)*SIND(THITA)*SIND(PHI)+ 
k  +COSD(PlrSIND(THITA)) 

P=SIND(Pl)*uOSD(THITA)*SIND(PHI)-COSD(Pl)*SIND(THITA) 

TEMP1=P*CEXP(J*DM)*CUR(G)/SIN(K0*LS(G)) 

TEMP2=CEXP(J*DM)*CUR(G)*SIND(P1)*COSD(PHI)/SIN(KO*LS(G)) 

IF(ABS(EM-K0)/K0.LT.l.E-2.OR.ABS(EM+K0)/K0.LT.l.E-2)THEN 

KG=LS(G)*SIN(K0*LS(G)) 

ELSE 

KG=-2*KO*(COS(EM*LS(G))-COS(KO*LS(G)))/(EM**2— 1\0**2) 

ENDIF 

NTHI0=TEMP1*KG+NTHI0 
NPHI0=NPHI0+KG*TEMP2 
58  CONTINUE 

CR=CABS(NTHI0)**2+CABS(NPHI0)**2 

KRON=(A*KY+B*KZ)/KX 

RC^K0V*2)*(UN**2)*CR/(4.W(ABS(A)**2+ABS(B)**2+ 
k  +ABS(kRON)**2)) 

IF((RC/((30./F)**2)).LT.l.E-6)THEN 
RCS  =  10.*ALOG10(l.E-6) 

ELSE 

RCS=10.*ALOG10(RC/((3G./F)**2)) 

ENDIF 

RETURN 

END 
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APPENDIX  B 
PROGRAM  TREE 


C  PROGRAM  TREE 

C 

C  PROGRAM  TO  GENERATE  A  PLANAR  FRACTAL  TREE  WITH 
C  REDUCTION  FACTOR  (COND)  AND  BRANCHING  ANGLE  THETA 
C  WRITTEN  BY  T.R.  NELSON,  PhD,  UNIVERSITY  OF  CALIFORNIA 
C  SAN  DIEGO,  LA  JOLLA,  CA,  90293,  AND  MODIFIED  BY 
C  LT.  JOHN  DEMIRIS  TO  GENERATE  THE  INPUT  DATA  FOR 
C  THE  RCS  PROGRAM. 

C 

C  INPUT  DATA: 

C 

C  L  =  NUMBER  OF  BRANCHING  LEVELS. 

C  N  =  NUMBER  OF  BRANCH  SEGMENTS. 

C  ILEN  =  INITIAL  LENGTH. 

C  ISP  =  INITIAL  STARTING  POINT. 

C  COND  =  REDUCTION  FACTOR. 

C  THETA  =  BRANCHING  ANGLE. 

C 

C  OUTPUT  DATA: 

C 

C  IX,  IY  =  COORDINATES  OF  FIRST  AND  END  POINT  OF  EACH 
C  GENERATED  BRANCH. 

C  INPUT  DATA  FOR  RCS  PROGRAM 

C 

REAL  7Z  (1:5,1: 1024), ISP, ILEN, IX, IY.ITH1 

OPEN  (UNIT=1,  FILE  =  'INGEOM'.FORM  =  'FORMATTED') 

OPEN  UNIT=2,  FILE  =  'OUT',  FORM  =  'FORMATTED') 

OPEN  (UNIT=10,FILE  ='INPUT1',F0RM='F0RMATTED') 
READ(1,*,END=10)  L,N, ILEN, ISP, COND, THETA 
10  CLOSE(l) 

ZZ(  1,1)  =  ISP 
ZZ  (4,1)  =  ILEN 
ZZ  (3,1)  =  0. 

ZL  =  ZZ  (4,1) 

ZZ  (2,1)  =  0. 

ZN  =  N 

PI  =  3.14159265 
IX  =  0. 

IY  =  0.-ILEN 
IX  =  ZZ(2,1) 

IY  =  ZZ(1,1) 

DO  100  LOOP  =1,L 
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c 

c 

c 

c 


(’ 

c 

c 

c 


c 

c 

c 

c 


ZLOOP  =  LOOP 
NPTS  =  N**ZLOOP 
INC  =  0 

NPREV  =  N**(ZLOOP-l) 

PREV  =  NPREV 
ANGLE  =  (ZN— 1.0)/2.0 
ZINDX  =  -ANGLE 
DO  200  J  =1,NPTS 
ZJ  =  J 

IF  (INC.GE.N)  ZINDX=-ANGLE 

IF  (INC.GE.N)  INC=0 

IW=NPTS-J+1 

IR  =  PREV— ZJ/ZN+ 1.0 

ZL=ZZ(4,IR) 

T  =  ZZ(3,IR) 

X  =  ZZ  (1  JR) 

Y  =  ZZ(2,IR) 

ZL1  =  ZL*COND 

IF  (ZL1.LT.0.1)  GOTO  300 

TH1  =  ZL1/33.0 

T1  =  T+THETA*ZINDX 

ZZ(3,1W)  =  T1 

T2  =  T1*PI/180.0 

IX  =  Y 

IY  =  X 

TEMPX1=IX 

TEMPY1=IY 


IX, IY  ARE  THE  COORDINATES  OF  THE  FIRST  POINT  OF  THE 
BRANCH  ALONG  THE  X.Y  AXIS  RESPECTIVELY 

WRITE(2,*)  IX, IY 
XI  =  ZLl*COS(T2)+X 
Y1  *  ZL1*SIN(T2)+Y 
IX=Y1 
IY=X1 


IX, IY  ARE  THE  COORDINATES  OF  THE  END  POINT  OF  THE 
BRANCH  ALONG  THE  X,Y  AXIS  RESPECTIVELY 


WRITE(2,*)  IX, IY 

TEMPX2-IX 

TEMPY2-1Y 


S  AND  T  ARE  THE  COORDINATES 
BRANCH 


OF  THE  CENTER  OF  EACH 


S  =  (TEMPX2+TEMPXl)/2. 

T  =  (TEMPY2+TEMPY1  )/2. 

PSI 1 = AT  A  N((TEM  PX2-TEM  PX  i  )/(TEM  PY2-TEM  PY 1 )) 
PSI-PSI1*180./P1 
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ZLHALF-ZL1/2. 

THlHALF=THl/2 

WRITEflO,*)  ZLHALF,TH1HALF,S,T,PSI 


ZZ( 

T,IW)=X1 

ZZ( 

2,IW)=Y1 

• 

ZZ( 

’4,IW)=ZL1 

zir* 

bX=ZINDX+1.0 

INC=1NC+1 

t 

200 

CONTINUE 

100 

CONTINUE 

300 

STOP 

END 
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